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SYNOPSIS 


In  the  First  Quarterly  Technical  Report,  14  May  1966,  a  functional  system 
description  of  a  LASA  Signal  Processing  System  was  hypothesized.  Certain 
system  parameters  were  outlined  and  bounded  to  permit  a  definition  of  the  signal 
processing  concept.  This,  the  Second  Quarterly  Technical  Report,  discusses 
important  system  parameters  in  greater  detail  to  provide  further  technical  sup¬ 
port  to  the  system  concept  and  describes  some  of  the  supporting  efforts. 

Sections  1  and  2  describe  the  work  accomplished  and  future  plans.  Appen¬ 
dices  A  to  C  contain  factual  detail  in  the  areas  of  beam  analysis,  signal  processing 
and  processing  systems.  Appendix  D  describes  effort  relating  to  communications 
instrumentation.  A  summary  description  of  the  computer  programs  currently 
under  development  is  given  in  Appendix  E.  Overall  program  comments  and 
recommendations  on  relevant  topics  of  seismological  interest  are  presented  in 
Appendix  F. 
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Section  1 


WORK  ACCOMPLISHED  AND  RESULTS  OBTAINED 

During  this  quarter,  work  supporting  the  overall  system  function  described 
in  the  First  Quarterly  Technical  Report*  has  been  accomplished.  The  scaling 
requirements  in  the  detection  processing  function  have  been  detailed;  process 
simulation  programs  have  been  prepared  and  limited  amounts  of  data  analyzed; 
gross  system  functions  both  internal  and  external  to  event  and  detection  proc¬ 
essing  have  been  partitioned  and  further  analyzed;  and  the  experimental  display 
has  been  constructed.  The  array  study  task  has  been  assigned  a  low  priority 
largely  because  of  the  desirability  of  testing,  by  means  of  the  simulation  programs 
and  the  display,  the  actual  criticality  of  seismometer  placement  and  the  real 
importance  of  side -lobe  reduction. 

To  facilitate  the  presentation  of  factual  detail  in  this  section,  the  following 
task  summaries  are  organized  according  to  the  order  delineated  in  the  contract. 

1.1  ARRAY  STUDY 

A  new  program  for  the  calculation  of  beam  patterns  is  in  preparation.  This 
program  will  facilitate  array  synthesis  by  making  use  of  superposition  of  array 
elements  and  corresponding  beam  patterns.  To  accomplish  this  will  require 
a  capability  to  plot  phase  as  well  as  logarithmic  amplitude  contours.  This  pro¬ 
gram  may  become  useful  for  possible  Montana  LASA  side-lobe  reduction  as  well 
as  further  study  of  regular  polygon  arrays. 

An  important  simulation  result  concerning  subarray  geometry  is  that,  under 
most  circumstances,  subarray  signal -to -noise  ratio  gain  can  be  enhanced  by 
not  using  the  B-ring  of  six  seismometers  of  the  present  subarrays,  thus, 


*IBM  First  Quarterly  Technical  Report,  "LASA  Signal  Processing,  Simulation 
and  Communications  Study,"  Contract  No.  AF  19(628)-5948,  May  27,  1966. 


decreasing  the  signal  processing  load.  This  result  appears  to  be  consistent  with 
results  of  others  who  have  found  that  it  would  be  advantageous  to  increase  the 
subarray  diameter.  It  has  also  been  found  that  one  of  the  six  subarray  beams 
may  be  eliminated  by  dispersing  the  subarray  centers.  Results  of  studies  per¬ 
taining  to  beam  requirements  are  presented  in  Appendix  A. 

1.2  STEERING  DELAY  STUDY 

Steering  Delay  Analysis  programs  have  been  modif  ied  and  are  ready  to  be 
applied  to  a  significant  number  of  LASA  wave-front  arrivals.  Central  to  these 
programs  is  the  Least  Squares  Wavefront  program  which  has  been  amended  so 
as  to  calculate  best  fitting  quadratic  wave  fronts  as  well  as  plane  fronts.  In  the 
event  that  the  second  order  (curvature)  parameters  should  contain  a  significant 
predictable  component,  this  addition  will  be  valuable  for  reducing  the  magnitude 
of  the  residuals.  Other  programs  which  are  to  be  used  in  the  Steering  Delay 
Study  are  Seismic  Steering  Delay  Anomalies,  Ray  Tracing,  Inverse  Velocity 
Space  Mapping,  and  Cross -Covariance. 

Limited  experience  to  date  in  measuring  delay  anomalies  as  well  as  experi¬ 
ence  from  analytical  and  processing  simulation  studies  indicates  that  steering 
delay  errors  amounting  to  a  beam  loss  of  a  few  db  at  1.5  Hz  may  be  inherent  in 
the  use  of  the  preformed  beams  for  detection  and  routine  event  processing  pur¬ 
poses.  For  this  reason,  a  special  processing  function  using  event -tailored  beams 
is  intended  for  those  detections  of  particular  interest.  Such  detections  will  also 
include  wave  arrivals  needed  for  updating  system  calibration. 

1.3  PROCESSING  REQUIREMENTS 

A  method  of  calculating  signal  scaling  in  the  beamforming-filtering  opera¬ 
tion,  so  as  to  retain  a  sufficiently  high  saturation  level  while  preventing  excessive 
deterioration  due  to  round-off  error,  has  been  developed  and  is  presented  in 
Appendix  B. 

A  limited  amount  of  simulation  of  subarray  beamforming  has  been  performed 
for  the  Longshot  event  with  several  types  of  filtering  (see  Appendix  B).  Thus,  a 
method  has  been  established  to  permit  identification  of  those  filters  which  achieve 
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maximum  signal -to -noise  ratio.  Furthermore,  an  anticipated  improvement  in 
subarray  signal -to -noise  ratio  obtained  by  dropping  the  B-ring  of  seismometers 
was  substantiated  by  the  simulation.  This  fact  will  lead  to  a  modest  saving  of 
computational  effort  required  to  implement  a  signal  processing  system  while 
generally  increasing  system  performance. 

It  has  been  determined  that  the  sampling  losses  associated  with  detection 
processing  are  acceptable  if  the  subarray  and  LASA  beamforming  are  executed 
at  a  rate  of  10  and  5  Hz,  respectively.  The  beamforming  phase  losses  are 
reasonably  contained  if  the  subarray  and  LASA  beamforming  resolutions  are 
20  and  10  Hz,  respectively.  In  this  manner,  the  subarray  beam  bandpass  filter¬ 
ing  may  be  performed  at  10  Hz  and  the  pre-detection  and  detection  processing  of 
LASA  beams  may  be  executed  at  5  Hz  (Appendix  B). 

Algorithm  specifications  have  been  refined  to  permit  more  efficient  micro¬ 
coding,  and  each  algorithm  is  being  simulated  in  assembly  language. 

A  working  breadboard  model  of  a  beam  display  has  been  built,  and  is  being 
used  to  evaluate  its  applicability  in  the  system  monitor  function  and  to  study 
beam  formations,  particularly  for  event  processing. 

Further  consideration  given  to  the  problems  of  event  processing  appear  to 
lead  in  the  direction  of  two  changes  in  emphasis.  First,  the  likelihood  that  pre¬ 
formed  beams  will  always  be  a  few  db  below  tailor-made  beams,  especially  at 
the  high  frequencies,  suggests  the  importance  of  a  capability  for  special  proc¬ 
essing  beyond  that  using  preformed  beams  routinely  generated  in  the  event  beam- 
former.  Secondly,  the  probable  existence  of  post-P  arrivals  and  side-lobe 
detections  for  larger  events  requires  a  capability  for  the  event  processor  to 
sort  and  identify  many  of  these  arrivals,  leading  to  both  logic  requirements  and 
additional  beamforming  capability  outside  the  P-teleseismic  zone. 

Considerable  work  has  been  done  in  defining  such  operational  system  functions 
as  monitor,  maintenance,  diagnostics,  interface  requirements,  operation  consoles, 
displays,  and  system  reliability.  Particular  emphasis  has  been  given  to  the 
repairable  duplex  system  mode  of  operation  to  assure  maximum  system  usability, 
especially  in  the  detection  processing  and  recording  functions.  The  techniques 
employed  in  these  processing  system  studies  are  presented  in  Appendix  C. 
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1.4  COMMUNICATIONS 


The  activity  in  communications  directed  at  the  transmission  of  LASA  data 
over  existing  common  carrier  systems  has  indicated  the  trade-offs  possible  in 
relating  cost  to  word  size,  sampling  rate  and  distance  between  the  array 
and  the  signal  processing  system.  A  description  of  this  work  is  included  for  the 
sake  of  completeness  in  Appendix  D. 

1.5  DATA  ANALYSIS  PROGRAMS 

A  number  of  programs  have  been  written  or  adapted  for  purposes  of  data 
analysis  and  simulation.  A  brief  description  of  their  function  and  input -output 
requirements  and  capabilities  is  presented  in  Appendix  E. 

1.6  COMMENTS 

The  IBM  LASA  program  has  been  reviewed  from  a  seismological  viewpoint 
by  the  project's  seismology  consultant.  A  summary  of  his  comments  is  pre¬ 
sented  in  Appendix  F.  He  has, on  the  whole, found  the  approaches  taken  in  devel¬ 
oping  the  functional  processing  system  to  be  reasonable  and  based  on  seismo¬ 
logical  facts  as  well  as  they  are  known  at  this  time.  He  has  also  expressed  an 
interest  and  encourages  that  further  seismic  research  be  carried  out  using 
LASA,  not  only  to  improve  the  system  performance  itself,  but  also  to  obtain  a 
more  detailed  understanding  of  the  earth's  interior. 
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Section  2 


PLANS 

The  following  tasks  are  planned  for  execution  during  the  remainder  of  the 
program.  The  tasks  are  discussed  in  approximate  order  of  priority.  Note  that 
the  ordering  will  be  dependent  upon  the  results  of  related  tasks  and  the  avail¬ 
ability  of  time  and  data. 

2.1  TIME  DELAY  STUDY 

A  revised  linear -quadratic  wavefront  program  will  be  used  to  process  a  few 
seismic  events  which  have  been  received  at  LASA  and  whose  arrival  times  at 
each  subarray  have  been  measured.  An  attempt  will  be  made  to  organize  these 
arrivals  on  a  world-wide  basis,  and  to  investigate  ways  to  make  the  steering 
delay  errors  reasonably  small. 

A  modest  program  of  measuring  arrival  times  from  LASA  event  tapes  will 
be  undertaken.  For  strong  events,  it  is  expected  that  one  of  several  possible 
methods  of  thresholding  will  be  acceptable.  For  weak  events,  the  question  exists 
whether  beam  gain  peaking  or  correlation  methods  might  provide  usable  arrival 
time  differences.  Hopefully,  a  peaking  program  can  be  written  which  answers 
this  question.  The  objective  is  to  develop  methods  of  automatic  arrival  time 
difference  measurement  to  determine  the  residual  steering  error  as  a  function 
of  measurement  technique  and  array  aperture. 

2.2  LOGIC,  EQUATIONS  AND  PROGRAMS  FOR  EVENT  PROCESSING 

A  logic  will  be  developed  to  sort  and  identify  the  manifold  arrivals  or 
detections  which  are  expected  to  be  made  by  the  detection  processing  function. 
These  are  expected  to  include  both  side-lobe  detections  and  subsequent  phases. 
The  routine  event  processing  function  will  include  this  logic  as  well  as  appro¬ 
priate  preformed  event  beams  and  associated  event  location  methods. 
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In  addition,  equations  will  be  formulated  and  programs  written  for  special 
processing  of  interesting  events.  These  programs  are  likely  to  include  adaptive 
beamformation  to  obtain  the  best  possible  steering  delays  for  selected  events. 

2.3  DISPLAYS 

A  breadboard  model  beam  display  has  been  constructed  and  is  ready  for 
experimental  use  on  taped  events.  Typical  experiments  which  can  now  be 
pursued  include:  display  of  subarray  "glitches”  as  differentiated  from  real 
events,  display  of  the  simulated  system  impulse  response,  effect  on  display  of 
various  combinations  of  time-delay  errors,  effect  of  various  display  intensity 
scales,  display  of  detection  as  well  as  event  beam  patterns,  operating  consider¬ 
ations  such  as  retention  of  detected  arrivals  until  they  are  turned  off,  evaluation 
of  various  slow-motion  time  scales  and  detection  techniques. 

2.4  DATA  ANALYSIS  ACTIVITY 

Results  of  subarray  data  analysis  performed  on  the  Longshot  event  are 
reported  in  Appendix  B.  This  analysis,  as  well  as  LASA  beam  analysis  and 
corresponding  filtering  and  its  effect  on  both  detection  and  event  processing 
performance,  and  other  seismic  events  analysis,  will  be  continued. 

2.5  MICROCODE  SIMULATION 

Simulation  of  the  algorithms  as  defined  in  Appendix  C,  with  a  view  to  attain¬ 
ing  maximum  speed  efficiency  in  microprogramming  and  simulation  of  possible 
implementation  techniques,  will  be  accomplished. 

2.6  DETECTION/EVENT  PROCESSING  CONTROL 

The  control  function  philosophy  for  a  duplex  signal  processing  system  will 
be  defined. 


2-2 


2.7  RAY  TRACING 


The  ray  tracing  program  may  be  modified  to  efficiently  calculate  post-P 
arrival  travel  times,  to  automatically  calculate  wave-front  curvatures  and 
intensity  for  comparison  with  measured  curvatures,  and  to  accept  velocity 
profiles  with  a  continuous  velocity  gradient  so  as  to  eliminate  any  spurious 
caustics. 

2.8  BEAM  PATTERNS 

A  beam  pattern  program  to  plot  amplitude  and  phase  contours  will  be 
investigated  for  possible  application  to  Montana  LASA  beam  pattern  improve¬ 
ment  and  regular  polygon  array  synthesis. 

2.9  LOCATION  CAPABILITY 

The  cross -covariance  program  can  be  applied  to  arrivals  at  world-wide 
seismic  stations  to  test  the  capability  of  providing  more  accurate  relative 
arrival  times  for  improving  the  location  accuracy  capability  of  a  potential 
world  network,  provided  available  records  exhibit  sufficient  timing  accuracy  to 
support  the  experiment. 
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Appendix  A 


BEAM  ANALYSES 


A.l  INTRODUCTION 

In  this  Appendix  the  activity  pertinent  to  beam  analysis  is  presented.  The 
program  described  in  Section  A. 2  has  been  employed  to  provide  the  velocity 
space  mapping,  and  its  description  is  included  herein  for  the  sake  of  complete¬ 
ness.  The  analysis  summary  in  Section  A. 3  describes  the  method  of  estimating 
beam-count  requirements.  The  dimensions  of  the  beams  are  discussed  in  Sec¬ 
tion  A. 4,  which  together  with  Section  A.5,  constitutes  a  baseline  estimate  for 
event  location  before  further  beam  processing.  Some  brief  remarks  concerning 
the  generation  of  beam  steering  delays  and  their  significance  with  respect  to  array 
gain  are  presented  in  Section  A.5.  Section  A. 6  contains  the  description  of  a  tech¬ 
nique  which  allows  a  reduction  in  subarray  beamforming  requirements  simply  by 
selective  steering  and  Section  A.7  concisely  addresses  location  estimates. 

A. 2  VELOCITY  SPACE  MAPPING 

The  program  used  for  preparing  inverse  velocity  space  maps  of  the  world 
land  ayeas  and  seismic  zones,  as  well  as  lines  of  constant  latitude  or  longitude, 
is  documented  in  this  section.  The  program  will  run  in  either  of  two  modes. 

In  Mode  1,  input  pairs  of  latitude  and  longitude  are  accepted  and  mapped  into 
inverse  velocity  coordinates.  In  Mode  2,  points  at  regular  intervals  of  latitude 
and  longitude  are  mapped  to  obtain  coordinate  lines.  World  map  latitude  and 
longitude  coordinates  used  for  preparing  the  inverse  velocity  map  were  supplied 
in  magnetic  tape  format  through  the  courtesy  of  the  IJ.  S.  Navy  Oceanographic 
Office  and  are  gratefully  acknowledged. 
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The  transformation  to  be  made  consists  of  calculating  the  range  and  azimuth 
by  means  of  spherical  trigonometry  and  then  converting  range  to  horizontal  phase 
velocity.  The  latter  transformation  is  based  on  the  ray  tracing  presented  in 
IBM's  LASA  Signal  Processing  Study.*  The  present  program's  intended  appli¬ 
cations  are  to  estimate  beam  counts,  to  roughly  predict  the  relationship  between 
geographic  locations  and  inverse  phase  velocity  (W),and  to  prepare  W-space  maps. 
It  does  not  contain  a  provision  to  correct  for  the  earth's  ellipticity. 

The  quadratic  formula  used  in  this  program  to  relate  range  to  velocity  is  a 
close  approximation  to  the  values  of  W  =  1  versus  XDR  derived  in  the  LASA 

Final  Report.* 

The  spherical  trigonometry  of  the  transformation  is  shown  in  Figure  A-l. 

The  range  angle  R  is  equal  to  side  c  of  the  spherical  tr  iangle  ABC,  and  can, 
therefore,  be  evaluated  from 

cos  R  =  sin  sin  q+  cos  eQ  cos  q  cos  (q>  -  <po),  (1) 


where: 


O  <  R  <  180°. 

The  azimuth  is  determined  by 

(2a) 

(2b) 


sin  AZI  = 


cos  q  sin(<p-<pQ) 
sin  R 


cos  AZI  = 


sin  e 


sin  6  cos  R 
o 


cosfl^  sin  R 


where: 

O  <  AZI  <  360°.  Both  equations  (2a  and 2b)  need  to  be  considered  to 
determine  the  correct  quadrant  of  AZI. 


*IBM  Final  Report,  "LASA  Signal  Processing  Study,"  Contract  No.  SD-296, 
July  15,  1965,  p.  B-9. 


A (<p,Q)  .  .  .  Point  at  Longitude  <p  (East),  Latitude  0  (North) 
(<p,  6  are  negative,  if  West  or  South) 

B(<Pq/  Qq)  .  .  .  Coordinates  of  Receiving  Array 
C  .  .  .  North  Pole 
R  ^  c  =  Range  angle 
AZI  =  Azimuth,  angle  CBA 


Figure  A-l.  Relation  of  Range  and  Azimuth  to 
Latitude  and  Longitude 


From  these  rules,  the  following  procedure  for  a  computer  program  to 
calculate  the  tan-1  function  with  output  between  -  7r  and  7ris  derived: 

a.  Range: 

1.  Calculate  cos  R  from  (1) 

2.  sin  R  =  'nJ  1  -  cos^  R 

3.  R.  =  tan-ijUBiL. 

hcos  R 

4.  If  cos  R  O,  R  =  R' 

5.  If  cos  R  <  O,  R  =  7r  -  R' 

b.  Azimuth 


1.  If  <p  =  (po : 

If  9  >  9  ,  AZI  =  O 
-  o 

If  9  <  9q,  AZI  =  7 r 

2.  If  (p  =  (PQ  ±  JT 

if  e  >  -  eQ,  azi  =  o 

If  6  <  -  0  ,  AZI  =  7r 
o 

3.  If  neither  <p=<p  nor  <p  =  cpQ  +  proceed  as  follows: 
If  e  =  f,  AZI  =  O 

if  e  =  -  azi  =  it 

If  0  ^  i  proceed  as  follows: 

If  sin  9  -  sin#  cosR  =  O,  AZI*  =  ^  . 

o  z 

If  sin  0  -  sin#Q  cosR  4  O,  proceed  as  follows: 


cos  9  cos  9  Q  sin  ((p-  (pQ) 

sin  9  -  sin  9  cos  R 
o 


A -4 


AZI'  =  tan 


If  sin  {(p 
If  sin  (<p 
If  sin  (cp 
If  sin  ((p 


(pQ)>  O  and  AZI'  >0,  AZI 
<pQ)  >  O  and  AZI'  <  O,  AZI 
<p  )  <  O  and  AZI*  >  O,  AZI 
<pQ)  <  O  and  AZI'  <  O,  AZI 


AZI' 

AZI'  +  n 
AZI'  +  7T 
AZI'  +  2tt 


A. 2.1  Range  and  Inverse  Horizontal  Phase  Velocity 

The  range  in  kilometers  is  the  range  angle  R,  multiplied  by  the  radius 
of  the  earth  Rg.  The  ray  tracing  to  which  reference  has  been  made  is  closely 
approximated  between  30°  and  95°  range  by  the  quadratic  formula 

RE  •  R  =  A  +  BW  +  CW2, 

where: 

(RE  •  R)  =  the  range  in  km 

W  =  the  inverse  phase  velocity  in  sec/km 
A  =  13,528 
B  =  -46,116 
C  =  -106. 

W  may  thus  be  calculated  (for  0.04  <  W  <  0.08)  from 

w  =  -w  [B 

A  plot  of  World  land  areas  with  longitude  and  latitude  coordinate  intersections 
is  shown  in  Figure  A-2. 
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Figure  A-2.  Inverse  Velocity  Space  Map  of  World  Land 

Areas  and  Seismic  Zones,  Centered  at  Billings, 
Montana 
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A. 3  BEAM  REQUIREMENTS 


To  further  check  Montana  LASA  beam  requirements  for  teleseismic  world 
coverage*  and  to  facilitate  beam  counts  for  any  other  possible  site,  a  set  of  film 
templates  of  close -packed  circles  has  been  prepared  as  shown  in  the  example  of 
Figure  A-3.  These  templates,  when  used  with  a  W-space  map  such  as  Figure  A-2, 
can  be  used  to  easily  determine  the  numbers  of  detection  beams  of  any  size 
required  to  cover  any  given  area  for  a  prescribed  array  site.  Templates  with 
circle  diameters  having  successive  ratios  of  ^TIT,  and  ranging  from  1/16  inch 
up  to  i  2  inches  have  been  made.  For  circles  smaller  than  1/16  inch,  the  appro¬ 
priate  numbers  can  be  estimated  by  dividing  the  map  area  in  W-space  by  the 
area  of  the  inscribed  hexagon  (for  nonduplicating  coverage)  belonging  to  a  single 
beam,  since  for  such  small  beams,  the  edge  effects  on  the  map  can  be  neglected. 

Figure  A-4  shows  beam  requirements  for  coverage  of  land  and  seismic 
areas,  as  well  as  land  only,  as  a  function  of  the  W-circle  diameter,  as  measured 
on  a  map  scaled  to  1  inch  =  0.02  sec/km. 

The  value  of  the  W  diameter  depends  upon  (1)  the  array  used  (e.g.,  21  sub¬ 
arrays  vs.  17  subarrays  vs.  13  subarrays),  (2)  the  maximum  loss  in  db  to  be 
tolerated,  and  (3)  the  frequency  at  which  this  loss  may  occur. 

To  obtain  the  beam  diameters  in  W  space  for  practical  beams,  reference  is 
made  to  beam  patterns  which  have  been  calculated  arid  presented  in  previous 
reports.  For  the  main  lobe,  which  is  of  concern  here,  the  following  quadratic 
approximations  for  W-space  beam  diameters  are  adequate: 

21  SA  LASA  beams:  D  =  0.0042  ''flj/f, 

17  SA  LASA  beams:  D  =  0.0070  ST/f, 

13  SA  LASA  beams:  D  =  0.0118  N~L/f, 

9  SA  LASA  beams:  D  =  0.0185  ^L/f, 

Subarray  beams:  D  =  0.105  NT/f, 


where: 


L  =  loss  in  db 

f  =  frequency  in  Hz 

D  =  beam  diameter  in  sec/km. 

♦IBM,  First  Quarterly  Technical  Report,  "LASA  Signal  Processing  Simulation 
and  Communications  Study,"  Contract  No.  AF  19(628)-5948,  May  27,  1966. 
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Figure  A-3.  Sample  Template  Showing  Close-Packed  Arrangement  of 
Beams  (Beam  Diameter,  D  =  0.5  Inch  -  Corresponding  to 
0.01  sec/km) 
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Number  of  Beams  to  Cover 


Beam  Diameter  in  W-space,  sec/km 


k  (km  1)^ 
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Figure  A-4.  Required  Number  of  Detection  Beams 
vs  Beam  Diameter  in  Inverse-Velocity 
Space  Centered  at  Billings,  Montana 
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As  a  sample  calculation,  let  it  be  required  to  determine  the  number  of  detec¬ 
tion  beams  needed  to  cover  land  and  seismic  areas,  and  permitted  to  lose  1.7  db 
(in  addition  to  a  1.3  db  subarray  beam  loss),  at  1.5  H:z,  using  all  21  subarrays. 

From  the  formula  for  21  subarrays,  the  W-space  beam  diameter  is 
D  =  0.0042  •  >IJ  1.7  /1.5  =  0.0036  sec/km.  From  Figure  A-4,  the  number  of 
such  beams  needed  to  cover  land  and  seismic  zones  from  Montana  is  N  =  900. 

This  number  was  previously  estimated  conservatively  at  1000.* 

A .4  BEAM  DIMENSIONS 

This  section  presents  an  estimate  of  the  accuracy  of  event  location  on  the 
earth's  surface,  in  terms  of  the  size  of  LASA  beams  in  inverse-velocity  (W) 
space. 

The  following  nomenclature  will  be  used  in  the  tables: 

=  inverse  phase  speed,  sec/km 
=  distance  in  range  direction,  km 
=  distance  in  azimuth  direction,  km 

km 

=  W  space  to  geography  scale  factors, 

=  from  ray  trace*  or 

=  b  +  2cW  from  quadratic  approximation 

Re  sin  9 
=  W 

=  diameter  of  beam  circle  in  W  space,  sec/km 

=  diameter  of  beam  circle  in  K  (wave  number)  space, 
cycles/km 

=  fDw,  where  f  =  frequency,  Hz 

*IBM,  First  Quarterly  Technical  Report,  "LASA  Sigral  Processing  Simulation 
and  Communications  Study,"  Contract  No.  AF  19(628)-5948,  May  27,  1966, 
p.  B-3. 
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R  =  range,  km 

2 

R  =  a  +  bW  ■+  cW  is  quadratic  approximation  to  range  vs. 
inverse  phase  speed  relation.  R  in  km 

0  =  range  angle,  radians  (or  0°,  in  degrees); 

0  =  R/R£  where  R^  =  radius  of  earth 

Dj^,  =  range  axis  and  azimuth  axis  of  beam  ellipse,  km 

dXR  dX. 

D  =  — —  D-D  =  — —  D 
R  dW  W’  A  dW  W 

D  =  larger  of  DR  or  D^ 
a  =  13,528;  b  =  -46,116;  c  =  -106;  E.£  =  6370 


Table  A-l  presents  the  scale  factors  for  several  ranges  relating  small  changes 
in  inverse  phase  speed,  W,  to  kilometers  on  the  earth's  surface,  in  the  range  and 
azimuth  directions.  A  small  single  frequency  constant  loss  circle  in  W-space 
will  transform  into  approximately  an  ellipse  in  real  coordinates.  The  scale 
factors  of  Table  A-l  are  the  ratios  of  the  range  and  azimuth  semi-axes,  respec¬ 
tively,  to  radius  of  the  W-space  circle.  The  scale  factors  of  Table  A-l  pre¬ 
sented  for  values  of  W  of  approximately  0.04,  0.06,  and  0.08,  have  been  obtained 
from  a  previously  published  ray  tracing.*  An  additional  value  of  W  =  0.0769  has 

dXR 

been  inserted  because  the  value  of  =  -409,800  at  W  =  0.0801  is  artificially 

high.  This  anomaly  appears  to  have  been  caused  by  a  "caustic"  generated  by  a 
first-derivative  discontinuity  in  the  piecewise  linear  velocity  profile  used  in  the 
ray  trace. 


*IBM  Final  Report,  "LASA  Signal  Processing  Study,"  Contract  No.  SD-296, 
July  15,  1965. 
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Table  A-l.  Range,  Range  Angle,  and  Scale  Factors  Obtained  From 
Ray  Trace,  Velocity  Profile  2 


w 

R 

e 

dXR 

dW 

dX. 

A 

dW 

0.0412 

10454 

94.0° 

177,200 

154,200 

0.0604 

6569 

59.1° 

131,600 

90,480 

0.0769 

4178 

37.6° 

191,300 

50,510 

0.0801 

3419 

30.8° 

409,800 

40,660 

In  Table  A -2,  the  same  quantities  are  presented  after  a  quadratic  fit  has 
been  made  to  the  Range  versus  W  relationship.  The  artificial  discontinuity  near 
W  =  0.08  has  thus  been  removed. 


Table  A-2.  Range,  Range  Angle,  and  Scale  Factors  From  Quadratic 
Approximation  to  Ray  Trace,  Velocity  Profile  2 


W 

R 

e 

dXR 

dW 

dxA 

dW 

0.04 

10083 

90.69° 

126,100 

159,200 

0.06 

7161 

64.41° 

131,600 

95,700 

0.08 

3439 

30.93° 

206,100 

40,900 

The  actual  velocity  profile  of  the  earth  may  have  significant  velocity 
discontinuities  at  larger  depths,  and  thus  caustics.  However,  presently  avail¬ 
able  data  is  insufficient  to  come  to  an  accurate  conclusion  on  this  question. 
Therefore,  for  the  present,  estimates  will  be  based  on  a  smooth  profile,  so  that 
typical  scale  factors  would  appear  to  range  up  to  about  200,000  km/(sec/km). 
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Table  A -3  presents  the  length  in  km  of  the  larger  double-semi -axis  of  the 
beam  ellipse  for  various  types  of  beams  and  the  three  values  of  W.  Three  db 
contours  of  1.5  Hz  detection  beams  are  evaluated  for  the  full  LASA  array 
(21  SA),  as  well  as  for  the  truncated  arrays  found  by  eliminating  the  outer  one 
or  two  rings  of  subarrays,  respectively.  Also  given  are  values  for  21  SA,  3  Hz, 
3db  event  beams  and  21  SA,  1.5  Hz,  3db  event  beams,,  Beam  diameters  meas¬ 
ured  in  this  manner  are  thus  seen  to  range  from  about  70  to  2100  km.  Detection 
beam  location  for  a  17  SA  beam  is  thus  of  the  order  of  300  to  1300  km  depending 
on  range,  whereas  event  beam  resolution  is  of  the  order  of  100  to  400  km. 

Table  A -3.  Diameters  for  Various  Beams,  Using  Scale  Factors  From 
Quadratic  Approximation  of  Table  A -2 


Type  of  Beam 

D  km 

SA 

L(db)* 

f 

°K 

Dw 

W  =  0.04 

W  =  0.06 

W  =  0.08 

21 

3.0 

1.5 

0.0073 

0.0049 

780 

470 

200 

17 

3.0 

1.5 

0.0120 

0.0080 

1270 

770 

330 

13 

3.0 

1.5 

0.0200 

0.0133 

2120 

1270 

540 

21 

3.0 

3.0 

0.0073 

0.0024 

380 

230 

100 

21 

1.5 

3.0 

0.0052 

0.0017 

270 

160 

70 

*L  (db)  is  the  loss  in  db  at  the  beam  circumference.  It  is  assumed  that  the 
full  loss  is  incurred  in  forming  the  LASA  beam,  i.e.,  the  subarray  beams 
are  formed  with  negligible  loss. 
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A. 5  ARRAY  STEERING  DELAYS 

The  loss  in  array  gain  in  db,  attributable  to  errors  in  steering  delays  is 
given  by 

2 

Loss  =  10  log  1Qe  (2irf a)  , 

where:  f  is  the  signal  frequency  in  Hz  and  cr  the  standard  deviation  of  steering 
errors  in  seconds,  i.e., 

I  <dV5tks>2’ 

k=l 

where:  6t^  =  the  relative  arrival  time  at  the  lc  subarray  center  and 

6t^g  =  the  steered  delay  at  the  k**1  subarray  center. 

At  a  detection  frequency  of  1.5  Hz,  a  standard  deviation  of  0.5  sec.  and  1.0  sec. 
in  the  steering  delays  would  produce  losses  of  1.0 db  and  3.8  db,  respectively. 
Therefore,  it  is  important  that  the  delays  of  even  the  presteered  detection  beams 
be  known  to  at  least  an  accuracy  of  0.1  sec,  with  a  standard  deviation  of  consid¬ 
erably  less  than  this  value.  To  get  the  best  possible  gain  at  event  beam  frequen¬ 
cies,  it  is  likely  that  some  adaptive  procedure  such  as  adjusting  the  arrival  time 
at  each  subarray  for  maximum  array  gain  (beam  peaking)  or  for  maximum  corre¬ 
lation,  will  be  required. 

The  method  by  which  the  detection  beam  delays  will  be  obtained  must  be 
essentially  an  empirical  one,  based  on  the  actual  arrival  times  of  actual  wave 
fronts  received  at  LASA.  It  is  expected  that  a  large  number  of  strong  event 
receptions  will  be  analyzed  by  means  of  a  Least  Squares  Wavefront  program  to 
determine  the  velocity  components,  curvatures,  and  the  deviations  of  actual 
arrival  times  from  plane  or  quadratic  wave  fronts.  It  is  hoped  that  when  the 
time  deviations  for  each  group  of  wave  fronts  are  averaged  and  compared,  the 
standard  deviations  will  be  sufficiently  small  to  permit  accurate  steering. 
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The  regions  on  the  earth  for  which  such  calibration  delays  can  be  measured 
will  generally  be  those  from  which  numbers  of  relatively  strong  seismic  events 
are  received.  For  other  regions,  from  which  only  weak  events  are  received,  an 
interpolation  method  will  be  required. 

Programs  for  analyzing  numbers  of  LASA  arrivals  are  now  being  completed 
and  the  task  of  analysis  will  be  initiated  shortly.  These  programs  will  calculate 
wave-front  directions,  velocities,  and  curvatures,  as  well  as  the  deviations  of 
actual  from  ideal  arrival  times.  The  program  assumes  the  existence  of  a  set  of 
arrival  times  for  each  subarray. 

The  arrival  time  inputs  to  the  Least  Squares  Wavefront  program  will 
initially  be  those  obtained  through  the  courtesy  of  the  Teledyne  Seismic  Data 
Laboratory.  We  intend  also  to  obtain  such  times  for  other  events  by  simple 
thresholding  on  strong  events,  and  by  cross -correlating  weaker  signals  to  obtain 
more  comprehensive  coverage. 

The  wave -front  deviations  which  have  been  obtained  by  the  above  procedure 
will  be  averaged  appropriately  to  obtain  sets  of  pre steered  beams  for  the  detec¬ 
tion  processor.  The  event  processor  is  expected  to  have  an  adaptive  capability 
to  adjust  the  steering  for  a  given  beam  so  as  to  maximize  array  gain  for  that 
event. 

A. 6  DISPERSION  OF  SURARRA.Y  BEAMS 

A.6.1  Discrete  Subarray  Beam  Steering 

In  subarray  beam  calculations  to  date,  it  has  been  assumed  that,  to  cover  a 
given  region  of  the  earth,  all  21  subarrays  are  aimed  at  the  center  of  the  region. 
The  resulting  subarray  beamforming  loss  is  then  zero  at  the  center,  and  1.3  db 
at  a  distance  of  6W  =  0.04  sec/km  from  the  center.  Each  of  the  21  subarrays 
has  its  own  set  of  seismometer  locations,  corresponding  to  its  rotation  relative 
to  subarray  El  or  C4,  each  of  which  has  a  long  arm  pointing  North,  Therefore, 
the  number  of  separate  subarray  beams  to  be  formed  is  21  N,  where  N  is  the 
number  of  subarray  beams  per  subarray,  or  the  number  of  regions  or  sectors 
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to  be  covered.  With  this  method  of  coincident  coverage,  N  has  to  be  equal  to  six 
for  world  coverage  at  no  more  than  1.3  db  loss.  Thus,  world  coverage  is  achieved 
by  using  126  separate  subarray  beams,  aimed  at  six  discrete  locations  and  losing 
from  0  to  1.3  db. 


A.6.2  Azimuth  Dispersion 

It  is  possible  to  smooth  out  this  loss  in  the  azimuth  direction  so  that  the  loss 
will  no  longer  be  as  little  as  zero  at  the  center,  but  in  return  it  will  be  no  greater 
than  0.64  db  anywhere  when  N  =  6.  The  corresponding  losses  for  N  =  5  or  4  are 
0.76  db  and  1.12  db,  respectively.  This  smoothing  is  achieved  by  dispersing  the 
21  N  subarray  beams  in  azimuth,  with  21  beams  uniformly  distributed  in  each 
sector  as  shown  in  Figure  A -5.  When  this  is  done,  the  maximum  subartay  beam¬ 
forming  loss  using  only  5  subarray  beamgfper  subarray  for  world  coverage  is 
0.5  db  less  than  when  6  subarray  beams  per  subarray  are  steered  discretely. 

If  P  is  a  point  at  range  R  to  be  covered  by  a  set  of  21  subarray  beams,  then 
the  21  subarray  beams  nearest  to  P  are  chosen.  This  will  include  one  beam 
for  each  subarray.  As  shown  in  Figure  A-5,  all  subarray  beams  are  centered 
at  a  range  c,  to  be  determined  for  optimum  coverage.  The  angle  2 <p  of  the 
sector  is  given  by 


2  (p  =  360°  (1) 

N 

where:  N  =  the  number  of  sectors  or  the  number  of  subarray  beams  per 
subarray. 

The  contribution  to  a  subarray  beam  amplitude  of  the  ith  subarray  beam  is 
SX.,  where  S  is  the  individual  signal  strength  and 


Here,  a  is  the  rate  of  subarray  beam  loss  fall -off  with  radius  in  inverse  velocity 

2 

space.  For  LASA  subarrays  at  a  1.5  Hz  detection  frequency,  a  =  93.5  (km/sec)  . 
As  shown  in  Figure  A-5,  ^  is  given  by 
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V 


Figure  A- 5.  Azimuth  Dispersed  Subarray  Beam 
Steering  in  Inverse -Velocity  Space 
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r. 

x 


(3) 


■r 


R2  +  c2  -  2Rc  cos  0 


0 


1/2 


The  beam  amplitude  is  thus  approximated  by  the  integral 


=  _L_  y^-a(R2+  c2-2Rc  cos0  )  d0 


(4) 


It  will  be  shown  that  for  practical  values  of  a,  R,  c  and  0,  the  exponent  is  suffi¬ 
ciently  small  that  it  may  be  replaced  by  the  constant  apd  linear  terms  in  its 
power  series  expansion.  Therefore,  we  may  write 


f> 


a  (R2  +  c2  -  2Rc  cos  0)1  d  0 


■] 


(5) 


=  1  -  a  (R2  +  c2  -  2Rc  -§y2’) . 


Moreover,  the  power  loss  resulting  from  subarray  beam  forming  is  thus 

L(R)  =  -20  1og1()A.  (6) 

The  value  of  subarray  beam  range  c  will  be  taken  such  that  the  loss  is  equalized 
at  the  extreme  ranges,  R  =  0.04  and  0,08  (sec/km),  or 

c  =  ibhp  *  \  <Rmin  +  Rmax>*  (?) 

The  greatest  loss  will  then  occur  at  R=  Rmin  and  R  =  Rmax,  and  the  losses  Lmax 
at  these  extreme  ranges  will  be  equal. 
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Table  A ^4  shows  the  value  of  <p  (Equation  1),  the  best  subarray  beam  range  c 
(Equation  7),  and  worst  loss  Lmax  (Equation  5,  6)  when  6,  5,  or  4  subarray  beams 
per  subarray  are  used  in  a  dispersed  mode  to  obtain  world  coverage;  a  =  93.5 
(km/sec)  for  1.5  db  detection  subarray  beams,  and  Rmin  and  Rmax  are  0.04  and 
0.08  sec/km,  respectively.  Column  5,  Lmin  will  be  discussed  in  Section  A.6.3. 

Table  A -4.  Numerical  Values  of  Best  Beam  Range  and  Worst  Loss 
for  6,5,  and  4  Subarray  Beams  Per  Subarray 


N 

(p  (rad.) 

c  (sec/km) 

Lmax  (db) 

Lmin  (db) 

6 

7t/6 

0.0628 

0.64 

0.28 

5 

tt/5 

0.0641 

0.76 

0.43 

4 

tt/4 

0.0666 

1.16 

0.70 

It  remains  to  show  that  the  linear  approximation  to  the  exponential  was  valid. 

2 

For  the  worst  case,  N  =  4,  a  =  93.5  (km/sec)  ,  the  largest  possible  value  of 
Ris  Rmax  =  0.0575  sec/km,  so  that  we  have 

-aR2 

e  maX  =  0.735,  whereas 


1-aR 


2 

max 


0.691 


This  worst  case  shows  that  the  linear  approximation  for  the  exponential  is  within 
about  6%  and  is  adequate  for  the  estimates  made  here. 


A.6.3  Possibilities  of  Further  Gains  by  Radial  Distribution 

To  determine  whether  additional  benefits  might  be  obtainable  by  further 
averaging  the  beam  losses  through  radial  distribution  of  the  subarray  beam 
centers,  we  first  look  whether  there  is  indeed  anything  to  average,  i.e.,  whether 
the  minimum  loss  along  a  radius  is  significantly  below  the  maximum  loss  at 
Rmin  and  Rmax.  As  shown  below,  a  very  marginal  improvement  might  be  possible. 
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The  minimum  loss  occurs  at  a  value  of  R  such  that 


^  =  -a  [2R  -  2c  ]  =  0,  or 

dR  cp 

R  =  c^  . 

(p 

At  this  value  of  R,  A  has  the  value 


(8) 


*  •  „  2  ,,  sin 

Amin  =  1  -  ac  (1 - g-3-] 

i,  and 

(9) 

<P 

Lmin  =  -20  log1Q  Amin. 

(10) 

Numerical  values  for  Lmin  are  shown  in  the  last  column  of  Table  A-4,  and 
are  in  the  order  of  0.3  db  better  than  Lmax.  Thus,  one  can  conclude  that  the 
addition  of  radial  dispersion  might  provide  a  further  reduction  of  perhaps  0.15  db 
in  the  subarray  beamforming  loss.  In  the  light  of  our  present  poor  knowledge  of 
some  much  more  important  factors  such  as  the  effects  of  timing  error,  this 
relatively  marginal  improvement  possibility  will  not  be  further  pursued  at  this 
time. 


A. 7  EVENT  LOCATION  ESTIMATE 

This  section  is  concerned  with  estimating  the  accuracy  with  which  a  seismic 
event  might  be  located  by  a  LASA  installation.  Consider  the  cases  of  routine 
processing,  and  special  processing  for  events  of  particular  interest.  An  event 
may  require  special  processing  because  its  travel  times  and  location  are 
desired  as  calibration  inputs  to  the  system,  or  because  of  scientific  or  other 
interest. 

In  routine  event  processing,  a  bundle  of  event  beams  is  centered  on  the 
appropriate  detection  beam  center  line,  and  used  to  locate  the  event.  In  this 
case,  the  following  sources  of  location  error  are  expected  to  accrue: 
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a.  Errors  due  to  finite  beamwidth 

b.  Errors  in  beamforming  due  to  incorrect  steering  delays 

c.  Errors  in  relating  measured  phase  speed  to  range. 

In  the  case  of  special  processing,  beamwidth  errors  (type  a)  will  not  exist  and 
an  attempt  is  made  to  minimize  steering  delay  errors  (type  b).  Errors  relating 
phase  speed  to  range  (type  c)  will  exist  for  both  routine  and  special  processing; 
their  magnitude  depends  entirely  upon  the  accuracy  to  which  the  source  location 
has  been  previously  calibrated. 

The  error  due  to  finite  beamwidth  (a)  is  closely  related  to  the  beam  resolution 
discussed  in  Section  A.4.  There,  estimates  of  beam  sizes  in  geographic  coordi¬ 
nates  were  obtained.  It  was  shown  that  event  beams  formed  from  21  subarrays 
and  losing  no  more  than  3  db  at  3  Hz  have  a  range  semi-axis  of  between  50  and 
190  km.  Thus,  one  can  expect,  with  routine  processing,  to  locate  an  event  with 
a  resolution  of  at  least  this  order,  provided  only  that  the  correct  beam  can  be 
selected  from  among  a  bundle  of  event  beams.  Because  the  main  part  of  the 
seismic  energy  is  contained  in  a  frequency  band  considerably  below  3  Hz,  and 
therefore,  a  number  of  beams  surrounding  the  best  one  will  generally  "light  up," 
two  methods  of  selecting  the  correct  inverse  velocity  space  components  suggest 
themselves. 

The  first  is  to  select  the  inverse  velocities  at  the  center  of  gravity  of  the 
bundle  of  event  beams.  This  method  is  computationally  simple  and  has  the 
advantage  of  retaining  the  low-frequency  beam  energy  while  at  the  same  time 
permitting  examination  of  the  side-lobe  structure.  Another  possible  method 
might  be  to  reduce  the  spread  of  the  event  beam  bundle  by  high-pass  filtering 
of  the  energized  beams.  The  disadvantage  of  this  would  seem  to  be  that  most 
of  the  beam  energy  which  is  at  lower  frequencies  would  be  disregarded. 

We  next  turn  our  attention  to  type  b  errors;  the  effect  of  incorrect  steering 
delays  leading  to  an  incorrect  value  of  the  inverse  phase  velocity  components. 

As  a  first  estimate  of  this  error, we  consider  just  two  seismometers,  separated 
by  an  array  aperture  of  200  km,  with  an  arrival  time  error  of  0.1  sec  committed 
at  one  of  them  for  a  broadside  wave  front.  At  an  inverse  phase  speed  of  0.04 
sec/km,  the  error  expressed  in  distance  would  be  0.1/0.04  =  2.5  km,  or  an  angle 
of  2.5/200  =  0.0125  radians  in  azimuth.  At  90°  range,  a  timing  error  of  0.1  sec 
would  thus  correspond  to  a  distance  of  approximately  80  km. 


A -21 


We  now  estimate  this  type  b  location  error  due  to  incorrect  steering  delays, 
when  using  a  LASA  array  rather  than  the  hypothetical  two  seismometers.  Assume 
N  subarrays  (N  =  21),  distributed  with  rotational  symmetry  and  with  distribution 
density  decreasing  linearly  with  distance  from  the  center,  out  to  a  radius 
R  (R  =  100  km).  We  further  assume  statistically  independent  arrival  time  errors 
with  standard  deviation  o  seconds.  Then,  the  RMS  distance  between  an  observed 
versus  the  real  source  location  in  Wx~W^  (inverse  velocity)  space  is  given  by 


AW 


=  0.00756  o  sec/km. 


The  actual  error  in  geographic  space  is  given  by 


AX  =  AW 


where  the  derivatives  are  those  given  in  Table  A -2,  Appendix  A.  The  value  of 

o 

the  derivative  is  shown  there  to  vary  in  absolute  value  from  40,900  km  /sec. 
to  206,700  km  /sec,  depending  on  range  and  depending  on  whether  the  range  or 
azimuth  directional  derivative  is  considered.  Thus,  the  distance  error  varies 
as:  AX=  a;cr  km,  where  the  values  of  ce  in  km/sec  are  given  as  a  function  of  W 
for  the  range  and  azimuth  directions  in  Table  A-5. 


Table  A-5.  Location  Error  Parameter  Values 


W  (sec/km) 

o;(km/sec) 

Range  Direction 

Azimuth  Direction 

0.04 

954 

1210 

0.06 

1000 

724 

0.08 

1560 

310 
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If  we  were  to  assume  an  absolutely  accurate  knowledge  of  the  arrival  times, 
and  consider  only  the  effect  of  the  quantized  sampling  interval  of  0.05  second, 
the  value  of  cr  would  be  1/  (40 ''(3")  seconds,  corresponding  to  AX  between  4.5 
and  27  km.  This  error  cannot  be  reduced  by  any  processing  scheme  save  a 
shorter  sampling  interval.  Assuming  arrival  time  error  standard  deviations  of 
c=  0.05  sec.  (optimistic)  and  a=  0.1  sec.  (pessimistic),  the  corresponding 
error  from  Table  A-5  would  thus  be  between  15  km  and  78  km,  and  between  31 
and  156  km,  respectively.  These  are  the  location  errors  which  might  reasonably 
be  expected  in  special  processing.  They  are,  as  expected,  less  than  the  errors 
incurred  from  routine  processing  by  preformed  event  beams. 

The  above  discussion  of  location  accuracy  for  a  single  LASA  is  now  briefly 
summarized  in  the  light  of  alternative  ways  to  process  events.  When  preformed 
event  beams  are  used  to  locate  events,  the  basic  timing  error  is  that  of  the  pre¬ 
formed  beams  compared  to  the  actual  wave  front  arriving,  and  the  resulting 
error  in  W-space  is  numerically  about  that  of  a  3  Hz  3db  beam  radius,  or 
perhaps  50  to  190  km.  In  this  case,  subsequent  beam  processing  such  as  cal¬ 
culating  centers  of  gravity  of  beam  clusters  should  help  materially  to  reduce 
this  source  of  error.  When  on  the  other  hand,  an  event  is  treated  individually 
with  emphasis  on  steering  a  beam  designed  adaptively  for  the  event,  then  the 
corresponding  least  squares  W-space  location  can  probably  not  be  still  further 
improved.  This  error  is  a  minimum  of  30  km  due  only  to  sampling  error,  but 
is  more  likely  to  be  50  km  to  100  km  for  actual  events.  To  both  methods  of 
processing  must  be  added  the  additional  error  of  W  to  geography  conversion, 
which  can  be  minimized  only  for  areas  where  calibration  events  have  been 
available. 
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Appendix  B 


SIGNAL  PROCESSING 


B.l  INTRODUCTION 

In  this  Appendix,  scaling  requirements  for  the  LASA  beamformer  are 
discussed. 

In  Section  B.2,  the  arithmetic  characteristics  of  the  beamformer  relevant 
to  the  topic  of  scaling  are  described. 

In  Section  B.3,  experimental  results  based  on  an  analysis  of  Longshot  data 
are  discussed.  These  results,  obtained  by  applying  various  bandpass  filters  to 
subarray  beam  outputs,  can  be  formulated  as  follows: 

a.  Omission  of  the  B-ring  in  each  subarray  will  enhance  the  conventional 
processing  by  about  1  db. 

b.  A  LASA  processing  gain  of  about  30  db,  relative  to  an  unprocessed 
seismometer  output,  can  be  anticipated  for  signals  with  dominant 
frequency  between  0.8  Hz  and  1.6  Hz. 

c.  The  seismic  noise  level  in  a  LASA  beam  can  be  as  low  as  0.03  nm. 

In  Section  B.4,  two  criteria  for  scaling  are  defined.  The  first  criterion 
demands  the  saturation  level  for  the  LASA  beamformer  to  be  _>  50  nm.  This 
ensures  that  all  saturating  signals  can  be  analyzed  with  ample  signal -to-noise 
ratio  by  means  of  the  "padded"  seismometers.  The  second  criterion  requires 
that  the  round-off  errors  (bias  as  well  as  random)  and  the  sampling  in  the  LASA 
beamformer  will  degrade  the  performance  by  an  acceptable  amount  (e.g.,  ldb 
or  less).  In  applying  the  second  criterion,  allowance  must  be  made  for  the  possi¬ 
bility  of  a  seismic  noise  level  in  a  LASA  beam  as  low  as  0.03  nm. 

In  Section  B.5,  scaling  and  sampling  for  two  beamformer  configurations  are 
examined.  To  satisfy  the  first  criterion,  a  total  accumulated  right-shift  of  5  bits 
must  be  performed  on  the  seismometer  outputs  prior  to  the  formation  of  LASA 
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beams.  Of  this  amount,  it  is  desirable  (and  possibly  necessary)  to  apply  a  1-bit 
right -shift  prior  to  the  formation  of  subarray  beams.  However,  the  latter  shift 
should  be  kept  as  small  as  possible  to  minimize  the  degrading  effects  of  the 
round-off  errors  on  the  LASA  beam  outputs. 

Two  configurations  are  considered.  In  one  configuration,  bandpass  filtering 
is  applied  to  the  subarray  beam  outputs  but  the  LASA  beams  are  not  filtered.  By 
compensating  for  the  (predictable)  bias  in  the  LASA  beams  generated  by  the 
round-off  errors  in  the  beamformer,  the  performance  loss  can  be  held  to  0.15  db 
(it  is  0.45  db  when  the  bias  is  not  compensated).  The  second  configuration  applies 
low-pass  filtering  to  the  subarray  beams  and  bandpass  filtering  to  the  LASA 
beams.  Here  the  performance  loss  due  to  round-off  errors  is  negligible  and 
compensation  for  round-off  bias  is  unnecessary. 

Both  configurations  utilize  a  rate  of  10  samples  per  second  for  the  formation 
and  filtering  of  subarray  beams  and  a  rate  of  five  samples  per  second  for  the 
formation  and  filtering  of  LASA  beams.  However,  the  subarray  beam  phasing  and 
the  LASA  beam  phasing  are  implemented  with  resolutions  corresponding  to  20 
and  10  samples  per  second,  respectively  (thus  limiting  the  phasing  loss  to  less 
than  0.35  db  for  signal  frequencies  up  to  1.6  Hz).  As  a  consequence,  the  post 
detection  sampling  rate  is  also  five  samples  per  second.  The  effect  of  this 
sampling  rate  on  the  beamformer  performance  can  be  reliably  determined  only 
by  experimental  means,  because  it  not  only  depends  on  the  length  of  the  post¬ 
detection  integration  interval,  but  also  on  the  detailed  characteristics  of  the 
signal  waveform.  Nonetheless,  a  rough  estimate  can  be  obtained  by  assuming 
the  signal  to  be  a  pure  sinusoid.  For  post -detection  integration  of  a  few  seconds, 
the  loss  in  performance  compared  to  analog  processing  will  be  less  than  0.5  db 
even  in  "worst-case"  situations  for  signal  frequencies  up  to  1.6  Hz. 

It  is  emphasized  that  the  beamformer  configurations  discussed  herein  are 
intended  primarily  as  typical  examples.  They  by  no  means  represent  final 
choices. 
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B.2  CONFIGURATION  AND  ARITHMETIC  CHARACTERISTICS 
OF  THE  LASA  BEAMFORMER 

Figure  B-l  provides  a  block  diagram  for  the  beamforming  process.  It  refers 
to  the  event  processor  as  well  as  the  detection  processor. 

Step  1  is  outside  the  beamformer  inasmuch  as  it  represents  the  process 
whereby  the  ground  motion  is  sensed  by  the  seismometer  and  is  converted  to  a 
digital  output  (word  length:  13  bits  plus  sign;  value  of  the  least  significant  bit: 

0.028  nm). 

The  beamforming  performed  in  steps  2  and  4  is  of  the  conventional  (delay  - 
and-sum)type.  The  sampling  rate  for  the  detection  processor  is  10  samples  per 
second  up  to  and  including  step  3.  It  is  proposed  to  drop  the  sampling  rate  to 
five  samples  per  second  thereafter,  provided  that  the  LASA  beam  phasing  be 
implemented  with  a  resolution  corresponding  to  10  samples  per  second.  As 
pointed  out  in  Section  B.5,  the  penalty  for  halving  the  sampling  rate  after  step  3 
is  acceptable  (on  the  order  of  0.5  db  or  less). 

For  each  step  in  Figure  B-l,  the  output  is  produced  in  digital  form.  Hence,  it 
will  possess  a  quantum  level, Q,  and  a  saturation  level,  SL.  Here  Q  is  defined  as  the 
physical  value  of  the  least  significant  bit  in  the  output,  while  SL  is  the  output  level 
above  which  overflow  may  have  occurred  in  the  arithmetic  operations  performed  in 
the  step.  Both  Q  and  SL  are  measured  in  nanometers  (nm)  and  are  always  referred  to 
the  seismometer  scaling  as  a  standard.  For  the  step  1  output:  Q  =  0.028  nanometers 
and  SL  =  (213-1)  Q  =  229  nm.  As  another  example,  consider  step  2,  the  subarray 
beamforming.  In  this  process,  a  number  of  seismometer  outputs  (e.g.,  K)  are 
added  to  form  the  step  2  output.  No  division  by  K  is  actually  performed.  Yet, 
we  shall  look  upon  this  step  2  output  as  the  average  of  the  K  seismometer  out¬ 
puts.  As  a  consequence,  the  quantum  level  for  the  step  2  output  is  Q/K,  where 
Q  denotes  the  step  1  output  quantum  level.  Because  the  beamformer  word 
length  of  15  bits  plus  sign  is  being  considered,  the  step  2  saturation  level,  SL, 
cannot  exceed  (2  -1)Q/K.  In  fact,  equality  holds  here  because  the  beamformer 

arithmetic  can  be  designed  to  keep  account  of  all  overflow  conditions  occurring 
during  the  addition  process  and  to  compensate  for  it  in  the  final  result.  Thus, 
effective  overflow  only  occurs  when  the  final  result  exceeds  the  register  length. 
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To  Threshold 


Figure  B-l.  Beamformer  Block  Diagram 
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We  emphasize,  that  in  this  example,  the  absence  of  shift  operations  in  step  2  has 
been  assumed.  Because  Q=0.028nm,  note  that  for  K  =  19,  the  step  2  quantum 
level  and  saturation  level  are  0.00137  and  48.3  nm,  respectively. 

As  noted  previously,  the  single  precision  word  length  in  the  beamformer  is 
15  bits  plus  sign.  When  "chopping"  a  double  precision  word  or  a  right-shifted 
single  precision  word, a  biased  round-off  error  is  incurred.  This  is  a  conse¬ 
quence  of  the  "two’s  complement"  arithmetic  used.  The  bias  is  always  negative 
and  has  its  magnitude  equal  to  half  of  the  least  significant  bit.  The  random  part 
of  the  round-off  error  is  taken  to  be  uniformly  distributed  between  plus  and 

minus  half  of  the  least  significant  bit.  Thus,  a  single  precision  word  with  quantum 

X 

level  Q,  when  shifted  to  the  right  by  X  bits,  will  possess  a  quantum  level  of  2  Q, 

X-i  X-1  i - 

a  bias  of  -2  Q  and  a  random  error  with  the  standard  deviation  2  Q/  N  3. 

The  effect  of  round-off  can  be  most  serious  for  recursive  filtering,  because 

these  operations  can  enormously  amplify  the  round-off  errors  generated  in  the 

filtering  process.  For  this  reason,  in  each  sampling  period,  all  products  of 

numbers  for  either  recursive  or  convolution  filtering  are  summed  using  their 

full  double  precision  word  length  and  only  the  final  result  is  chopped  to  single 

precision.  As  a  consequence,  at  each  sampling  instant,  only  one  round-off  error 

is  generated  by  the  filtering  operation. 

Denoting  the  quantum  level  of  the  input  to  a  recursive  filter  by  Q,  the  filter 

output  will  have  a  bias  of  the  form  BQ/2  and  a  standard  deviation  of  the  form 

Q  Nd/12.  For  a  filter  of  fixed  degree,  the  values  of  B  and  D  tend  to  increase 

rapidly  with  the  sampling  rate  f  (the  degree  of  a  recursive  filter  is  the  degree 

-1  s 

of  the  rational  function  in  Z  describing  the  filter).  Table  B-l  lists  the 
values  of  B  and  D  for  a  Butterworth  filter  of  degree  six  and  with  a  pass  band 
from  0.7  Hz  to  1.4  Hz.  Note  that  the  amplification  of  the  internally  generated 
round-off  errors  for  a  bandpass  filter  can  sometimes  be  substantially  suppressed 
by  implementing  the  filter  as  a  low-pass  filter  running  at  a  high  sampling  rate, 
followed  by  a  high-pass  filter  running  at  a  low  sampling  rate. 
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Table  B-l.  Filter  Amplification  Factors 


fs(Hz) 

A 

B 

C 

D 

5 

0 

O 

• 

oo 

! 

0.28 

5 

10 

0 

28.0 

0.14 

2,020 

20 

0 

1394.0 

0.07 

2,300,000 

A  recursive  filter  also  propagates  and  amplifies  bias  and  random  error 
already  present  in  the  filter  input.  The  corresponding  amplification  factors 
will  be  labeled  A  and  C  (see  Table  B-l).  If  m  and  cr  denote  input  bias  and  stan¬ 
dard  deviation  of  the  input  random  error,  then  Am  and  a  are  the  correspond¬ 
ing  quantities  for  the  propagated  error  in  the  filter  output.  In  Table  B-l,  C  has 
been  calculated  on  the  assumption  that  the  random  input  errors  are  uncorrelated 
from  one  sampling  instant  to  another.  When  this  condition  is  not  satisfied,  an 
upperbound  for  C  is  furnished  by  the  maximum  value  of  | H(f )  |^,  where  H(f) 
denotes  the  complex  (amplitude -phase)  frequency  response  of  the  filter.  For 
the  filters  to  be  considered  for  steps  3  and  5  in  Figure  B-l,  this  maximum  value 
is  unity,  hence,  C  <  1. 

B.3  EXPERIMENTAL  RESULTS 

To  obtain  some  insight  into  the  effectiveness  of  filtering  as  a  means  for 
increasing  the  signal -to -noise  ratio  of  LASA  beam  outputs,  a  short  computa¬ 
tional  analysis  has  been  made  of  some  of  the  Longshot  data.  The  computations 
were  performed  in  floating-point  arithmetic  using  an  IBM  7094. 

The  25  recorded  seismometer  outputs  of  subarray  B1  were  delayed  and 
summed.  The  sum  was  divided  by  25  to  obtain  the  "complete"  subarray  beam 
output.  In  addition,  the  19  seismometer  outputs,  remaining  after  omitting  the 
six  seismometers  in  the  B-ring  of  subarray  Bl,  were  likewise  delayed  and 
added.  Their  sum  was  divided  by  19  to  provide  a  "partial"  subarray  beam.  The 
phasing  of  the  seismometer  traces  prior  to  their  addition  was  implemented  with 
a  resolution  corresponding  to  20  samples  per  second.  However,  both  beam  out¬ 
puts  were  calculated  corresponding  to  a  real-time  rate  of  10  samples  per  second. 
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The  complete  beam,  the  partial  beam  and  the  output  of  the  center  seismom¬ 
eter,  Aq, of  subarray  B1  were  passed  through  the  same  digital  recursive  filter. 
The  three  outputs  so  obtained  were  used  to  calculate  signal-to-noise  ratios  in 
the  following  manner.  The  noise  power,  N,  was  calculated  by  averaging  the 
squares  of  850  successive  samples  preceding  the  signal  onset  for  each  of  the 
three  channels.  Thus,  N  is  the  average  noise  power  over  an  85-second  time 
period  immediately  preceding  the  signal.  The  average  signal  power,  S,  for 
each  of  the  three  channels  was  calculated  over  an  interval  of  2.6  sec.  This 
interval  was  chosen  to  yield  maximum  average  signal  power, S.  The  signal-to- 
noise  ratio  of  each  channel  was  defined  as  101og10  (S/N).  Note  that  this 
definition  ignores  the  post-detection  integration  gain. 

The  filter  gain  for  the  seismometer  was  defined  relative  to  the  unfiltered 
A()  seismometer  trace.  The  array  gain  for  a  filtered  "complete"  beam  or 
"partial"  beam  was  computed  relative  to  the  A^  trace  after  passing  it  through 
the  same  filter  as  that  used  for  the  beam.  The  total  gain  was  taken  as  the  sum 
of  filter  and  array  gains.  The  seismic  noise  standard  deviation  was  defined 


The  results  of  these  calculations  can  be  found  in  Table  B-2  for  a  number  of 
bandpass  filters.  The  filters  are  listed  by  their  3-db  pass  band,  preceded  by 
the  designation  "Bu"  (for  Butterworth)  or  "Ch"  (for  Chebychev).  All  filters 
were  of  degree  six.  The  Chebychev  filters  were  implemented  with  a  0.5 -db 
ripple. 

Note  that  the  "partial"  beam  outperforms  the  "complete"  beam  by  about 
1  db,  except  for  the  Chebychev  filter  with  pass  band  from  0.9  Hz  to  3.0  Hz. 

The  explanation  lies  in  the  fact  that  the  low-frequency  microseismic  noise  is 
strongly  correlated  among  the  seven  seismometers  located  at  the  center  and  on 
the  B-ring  of  subarray  Bl.  Hence,  the  array  gain  for  the  complete  beam  tends 
to  be  inferior  to  that  of  the  partial  beam  for  those  filters  which  pass  more  low- 
frequency  than  high-frequency  noise.  Thus,  it  is  recommended  that  the  subarray 
beams  be  formed  from  all  subarray  seismometers  except  those  in  the  B-ring. 
This  procedure  will  result  both  in  enhanced  signal-to-noise  ratio,  as  well  as  in 
a  lighter  computational  load. 
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Table  B-2.  Experimental  "Longshot"  Subarray  B1  Data 


Filter  Type 
and 

Pass  Band 
(Hz) 

Noise  Standard  Deviation  (nm) 

Filter 
Gain  (db) 
for  A0 

Array  Gain  (db) 

Total  Gain  (db) 

AO 

Seismometer 

Complete 

Beam 

Partial 

Beam 

Complete 

Beam 

Partial 

Beam 

Complete 

Beam 

Partial 

Beam 

Unfiltered 

2.77 

2.04 

1.77 

0.0 

2.0 

2.9 

2.0 

2.9 

Bu  0.6-2. 0 

0.89 

0.42 

0.38 

9.3 

5.9 

6.6 

15.2 

15.9 

Ch  0. 9-3.0 

0.67 

0.26 

0.27 

11.4 

7.5 

7.1 

18.9 

18.5 

Bu  0. 7-1.4 

0.73 

0.29 

0.24 

10.7 

7.4 

8.6 

18.1 

19.3 

Bu  0.8 -1.6 

0.64 

0.24 

0.21 

12.1 

7.9 

8.9 

20.0 

21.0 

Bu  0.8-1.4 

0.61 

0.22 

0.19 

12.4 

8.0 

9.2 

20.4 

21.6 

Bu  0.9-1.6 

0.51 

0.19 

0.17 

13.9 

8.0 

8.7 

21.9 

22.6 

Bu  0. 9-1.4 

0.46 

0.17 

0.15 

14.6 

7.9 

9.1 

22.5 

23.7 

Ch  0.9 -1.4 

0.44 

0.17 

0.14 

14.7 

7.9 

9.2 

22.6 

23.9 

As  is  to  be  expected,  Table  B-2  demonstrates  that  the  signal-to-noise  ratio 
tends  to  increase  as  the  filter  bandwidth  is  made  narrower.  However,  it  must 
be  remembered  that  different  signals  will  possess  different  bandwidths.  Thus, 
as  the  filter  bandwidth  is  decreased  to  match  the  signal  bandwidth,  the  possi¬ 
bility  exists  that  the  number  of  filters  per  LASA  beam  required  to  accommodate 
the  received  signals  may  have  to  be  increased  above  unity.  Hence,  a  compromise 
must  be  struck  between  the  signal-to-noise  ratio  to  be  achieved  and  the  total 
number  of  filters  to  be  implemented.  For  that  purpose,  a  preliminary  analysis 
was  made  of  data  supplied  by  SDI?  and  listing  location,  magnitude  and  dominant 
signal  frequency  for  more  than  100  events  received  by  the  Montana  LASA.  All 
dominant  frequencies  were  within  the  band  from  0.6  Hz  to  2.0  Hz.  More  than 
85  percent  were  in  the  band  from  0.8  Hz  to  1.6  Hz. 

According  to  Table  B-3,  the  Butterworth  filter  with  a  bandwidth  from  0.8  Hz 
to  1.6  Hz  outperforms  the  Butterworth  filter  with  a  bandwidth  from  0.6  Hz  to 
2.0  Hz  by  about  5  db.  It  is,  therefore,  tempting  to  conjecture  that  the  bandwidth 
from  0.8  Hz  to  1.6  Hz  will,  in  fact,  be  adequate  for  the  detection  processor. 

For  example,  one  may  speculate  that  decreasing  the  lower  corner  frequency 
may  not  be  necessary  for  two  reaons:  (1)  there  is  a  tendency  for  events  with  low 
dominant  frequency  to  be  also  strong  events;  (2)  because  the  microseismic  noise 
level  increases  rapidly  below  0.8  Hz,  decreasing  the  lower  corner  frequency 
may  be  ineffective  as  a  means  of  increasing  the  LASA  beam  signal-to-noise  ratio 
for  events  with  low  dominant  frequencies.  At  any  rate,  to  decide  these  issues,  a 
computational  analysis,  similar  to  that  reported  for  Longshot,  should  be  applied 
to  recorded  events  with  low  dominant  frequencies  as  well  as  to  events  with  high 
dominant  frequencies. 

Pending  such  an  analysis,  the  possibility  must  be  faced  that  the  use  of  a  single 
bandpass  filter  will  provide  inadequate  signal-to-noise  ratio  for  LASA  beam 
outputs.  For  this  reason,  two  specific  configurations  are  described  and  analyzed 
in  Section  B.5  for  the  purpose  of  determining  a  suitable  scaling  and  assessing 
the  loss  in  performance  due  to  round-off  errors. 

The  first  configuration  applies  filtering  to  subarray  beams  only.  The  LASA 
beams  are  not  filtered,  i.e.,  the  step  5  filtering  (see  Figure  B-l)  is  omitted. 
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The  step  3  filter,  running  at  10  samples  per  second,  must  reject  both  high- 
frequency  and  low-frequency  noise.  Thus,  a  bandpass  filter  would  seem  to  be 
a  suitable  choice. 

In  the  second  configuration,  filtering  is  applied  to  subarray  beams  (step  3) 
as  well  as  to  LASA  beams  (step  5);  the  latter  filtering  running  at  five  samples 
per  second.  For  example,  the  step  3  filtering  could  be  (but  need  not  be)  of  the 
low-pass  type,  designed  to  prevent  "frequency  folding"  of  the  seismic  noise 
when  the  sampling  rate  is  subsequently  dropped  to  five  samples  per  second. 

The  step  5  filtering  could  be  of  the  bandpass  type,  allowing  LASA  beams  to  be 
filtered  according  to  their  individual  needs  (if  required,  using  more  than  one 
filter  per  LASA  beam). 

The  computational  analysis  of  Longshot,  reported  on  above,  has  been  incom¬ 
plete  in  the  sense  that  the  LASA  beam  output  signal -to -noise  ratio  has  not  been 
experimentally  determined.  Experimental  results  obtained  by  the  Lincoln 
Laboratories  of  MIT  indicate  that  the  additional  array  gain  resulting  from  the 
noise  suppression  only,  can  be  calculated  on  the  conventional  assumption  of 
statistical  independence  of  the  noise  at  different  subarrays,  thus  providing  an 
added  array  gain  of  10  log^M  dbover  and  above  the  array  gain  of  subarray 
beams.  Here,  M  is  the  number  of  subarrays  used  for  the  LASA  beamforming. 

From  this  must  be  subtracted  the  signal  loss  due  to  incoherence  (probably  on 
the  order  of  1  db),  as  well  as  the  more  severe  signal  loss  due  to  the  usage  of 
incorrect  steering  delays.  This  second  kind  of  loss  can  occur  for  LASA  beams 
steered  to  regions  of  low  seismic  activity.  It  tends  to  be  proportional  to  the  square 
of  the  dominant  frequency  and  could  conceivably  run  as  high  as  2  db  or  more.  Allow¬ 
ing  a  total  of  3  db  for  these  signal  losses  and  taking  M  =  17,  it  follows  that  the  total 
gain  for  a  LASA  beam  can  be  as  high  as  21  +  (10  logjQl7)-3  =  30  db.  (Here,  21  db 
is  assumed  to  be  delivered  by  the  subarray  gain.  See  Table  B-2.) 

As  pointed  out  earlier,  the  post -detection  integration  gain  has  been  ignored 
in  the  preceding  considerations.  For  the  computational  analysis  of  Longshot, 
reported  on  above,  a  post -detection  interval  of  2.6  seconds  was  used  to  compute 
the  maximum  average  signal  power,  S.  This  somewhat  arbitrary  choice  of  post¬ 
detection  integration  time  can  be  improved  upon  by  noting  that,  for  maximum 
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signal-to-noise  ratio,  the  post -detection  integration  time,  T,  should  be  selected 
to  maximize  the  expression  '^"t"  S  (this  assertion  assumes  the  seismic  noise 
to  be  a  stationary  process  with  a  correlation  time  which  is  small  compared  with 
T).  It  is  probably  worthwhile  to  plot  S  versus  T  for  a  number  of  recorded 
seismic  events  for  the  purpose  of  selecting  a  value  of  T  suitable  for  all  seismic 
events,  i.e.,  entailing  an  acceptable  loss  in  signal -to-noise  ratio  for  the  majority 
of  signals.  In  the  situation  where  such  a  compromise  value  cannot  be  established, 
it  may  be  necessary  to  use  several  post -detection  integration  times  in  the 
detection  processor. 

B.4  SCALING  CRITERIA 

It  is  the  purpose  of  this  section  to  examine  the  scaling  for  the  beamforming 
process  as  it  has  been  defined  in  Section  B.2.  Such  a  scaling  must  satisfy  the 
following  criteria. 

a.  Criterion  (1)— The  saturation  level  of  the  beamforming  process  should 
be  high  enough  to  allow  adequate  detection,  localization  and  processing  of  the 
saturating  signals  by  means  of  the  "padded"  seismometers  only  (one  such 
seismometer  per  subarray). 

b.  Criterion  (2)— The  round-off  errors  generated  by  the  beamforming 
process,  as  well  as  the  sampling  rates  employed,  should  degrade  the  LASA 
beam  output  signal-to-noise  ratio  by  only  a  fraction  of  one  db. 

To  satisfy  criterion  (1)  we  shall  require  a  saturation  level  for  the  beam- 
former  of  about  50  nm  or  more.  Because  the  unfiltered  noise  level  at  the  LASA 
site  is  about  10  nm  or  less,  the  signal-to-noise  ratio  of  an  unfiltered  saturating 
event,  as  measured  by  a  "padded"  seismometer,  will  be  on  the  order  of  14  db  or 
more.  This  can  be  increased  to  a  signal-to-noise  ratio  of  more  than  20  db  by 
filtering  using  a  wide  pass  band  (e.g.,  from  0.6  Hz  to  2.0  Hz,  see  Table  B-2). 
Such  a  signal-to-noise  ratio  permits  detection  of  saturating  signals  on  a  single 
"padded"  seismometer  and  affords  the  feasibility  of  accurate  timing  for  event 
localization.  Hence,  a  saturation  level  of  about  50  nm  or  more  is  indeed 
satisfactory. 
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To  obtain  a  useful  quantitative  formulation  for  criterion  (2),  note  that  the 
bias,  m,  and  the  random  error,  e(t),  generated  in  a  LASA  beam  as  a  result  of 
round-off  errors  in  the  LASA  beamformer,  will  cause  a  loss  in  signal -to-noise 
ratio  of  the  LASA  beam  output.  This  loss  is  given  by  the  following  expression, 


Loss  =  5  log 


10 


r 


(<ftn  +  ^e)2df  +  2m2  Ki(0)  +  ^e(0)) 


I 


6o 


2  df 
n 


db.  (1) 


Here  $  and  <p  are  the  power  spectra,  respectively,  of  e(t)  and  of  the  seismic 
noise  n(t),  in  the  LASA  beam  output.  The  derivation  of  equation  (1)  is  based  on 
the  following  assumptions:  the  total  noise,  n(t)  +  e(t),  is  a  Gaussian  stationary 
process  with  zero  mean  and  with  a  correlation  time  which  is  short  compared  to 
the  duration  of  the  post -detection  interval. 

From  equation  (1),  it  follows  that 


<T  ^  / 

Loss  «  4.3  -^2  ll  + 


(2) 


where  a  and  cr  are  the  standard  deviations  of  e(t)  and  n(t),  respectively.  The 
approximation  (equatioir2)  is  based  on  the  following  assumptions:  (1)  the  spectra 
of  n(t)  and  e(t)  are  flat  over  their  respective  effective  frequency  bands;  (2)  the 
frequency  band  for  e(t)  extends  to  zero  frequency,  while  that  of  n(t)  does  not, 
i.e.,  $n(0)  =  0;  (3)  the  frequency  band  of  n(t)  is  included  in  that  of  e(t);  (4)  none¬ 
theless,  the  ratio  of  these  two  frequency  bands  can  be  approximated  by  unity; 

(5)  <jJcrn  is  appreciably  less  than  unity.  These  assumptions  are  rather  drastic 
in  view  of  the  real  situation.  They  tend  to  make  the  approximation  (equation  2) 
conservative,  i.e.,  usage  of  equation  2  tends  to  overestimate  the  loss  caused  by 
the  round-off  errors  in  the  beamformer. 
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Henceforth,  equation  (2)  will  be  applied  when  examining  whether  the  scaling 
of  the  beamforming  process  indeed  complies  with  criterion  (2).  For  these  appli¬ 
cations,  it  is  necessary  to  have  an  estimate  for  a  .  According  to  Table  B-3, 
the  noise  standard  deviation  for  a  subarray  beam  can  be  as  low  as  0.14  nm. 

Hence,  for  a  LASA  beam,  using  all  21  subarrays,  the  noise  can  be  as  low  as 
0.14/  ^21  =  0.03 nm.  Henceforth,  when  applying  equation  (2),  we  shall  assume 
c  =  0.03  nm. 

B.5  SCALING  CONSIDERATIONS 

In  this  section,  the  scaling  requirements  and  the  effects  of  round-off  errors 
are  discussed  for  two  configurations  of  the  LASA  beamformer.  Referring  to 
Figure  B-l,  in  the  first  configuration,  the  step  3  filter  is  a  sixth  degree  Butter- 
worth  filter  with  its  3-db  pass  band  from  0.8  Hz  to  1.6  Hz.  The  step  5  filter  is 
omitted,  i.e.,  the  step  4  output  and  the  step  5  output  are  identical.  In  the  second 
configuration,  the  step  3  filter  is  a  low-pass  third  degree  Butterworth  filter  with 
its  3-db  cut-off  at  2  Hz.  The  step  5  filter  is  taken  to  be  a  sixth  degree  Butter- 
worth  filter  with  the  3-db  pass  band  from  0.7  Hz  to  1.4  Hz.  The  step  3  filter 
for  both  configurations  is  implemented  to  run  at  a  rate  of  10  samples  per  second. 
The  step  5  filter  is  intended  to  run  at  5  samples  per  second. 

As  already  explained  in  Section  B.3,  the  first  configuration  can  be  used  if 
one  single  filter  applied  to  all  LASA  beams  is  capable  of  providing  adequate 
detection  capability.  In  the  eventuality  that  such  is  not  the  case,  and  that  differ¬ 
ent  LASA  beams  require  different  filters  or  that  a  LASA  beam  must  be  monitored 
using  more  than  one  filter,  the  first  configuration  cannot  be  employed  and  an 
implementation  similar  to  the  second  configuration  may  be  necessary. 

We  begin  by  discussing  the  requirement  imposed  by  criterion  (1)  formulated 
in  Section  B.4.  Disregarding  for  the  moment  the  filtering  performed  at  step  3 
(see  Figure  B-l)  and  assuming  that  M  subarrays  are  used  to  form  a  LASA  beam, 
it  follows  that  the  step  4  quantum  level  is  2  Q/KM  and  that  its  saturation  level 
is  (2  '  -1)2  Q/KM.  Here.Q  is  the  step  1  quantum  level,  i.e.,  Q  =  0.028  nm,  and 
K  is  the  number  of  seismometers  used  per  subarray.  The  integer  X  is  the 
total  number  of  bits  shifted  from  the  step  1  output  up  to  the  point  where 
the  delay -and-sum  process  is  performed  in  step  4  to  produce  the  LASA  beam.  A 
positive  value  of  X  indicates  a  right  shift. 
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We  assume  that  M  is  chosen  between  17  and  21,  and  that  K  lies  between  19 

(when  the  B-ring  is  omitted,  see  Section  B. 3)  and  25.  Hence,  the  saturation  level 

X  X 

for  step  4  will  lie  between  1.75  x  2  and  2.84  x  2  nm.  ForX  =  5,  these  values 
are  56  and91nm.  Criterion  (1)  in  Section  B.4  demands  a  saturation  level  of 
about  50  nm  or  more.  Hence,  a  total  accumulated  right-shift  of  5  bits  is  required 
to  satisfy  criterion  (1)  for  the  step  4  processing. 

This  total  right-shift  of  5  bits  only  ensures  that  criterion  (1)  is  met  for  the 

step  4  processing.  However,  it  also  should  be  satisfied  for  all  other  steps  as 

15  v 

well.  Specifically,  for  the  step  2  output,  the  saturation  level  is  (2  -1)2JQ/K  nm, 

where  y  is  the  number  of  bits  of  right-shift  applied  to  the  step  2  input  prior  to 
the  delay -and-sum  process.  For  K  =  19,  this  saturation  level  is  49  x  2^nm  and 
for  K  =  25,  it  is  37  x  2ynm.  Hence,  to  satisfy  criterion  (1),  no  right -shift  is 
needed  for  the  step  2  input  when  K  =  19,  while  a  right-shift  of  one  bit  is  adequate 
when  K  =  25.  Of  course,  for  larger  right -shifts  criterion  (1)  will  still  be  satis¬ 
fied.  However,  it  will  be  shown  below  that  minimization  of  the  effect  of  round¬ 
off  in  the  LASA  beam  output  (i.e.,  in  the  step  5  output)  demands  a  minimum  value 
of  y,  i.e.,  demands  a  minimum  right-shift  of  the  step  2  input. 

This  brings  us  to  the  filtering  in  steps  3  and  5  (see  Figure  B-l).  Note  that 
the  filters  to  be  implemented  for  the  two  configurations  all  have  a  complex 
(amplitude-phase)  frequency  response  whose  absolute  value  nowhere  exceeds 
unity.  This  fact,  together  with  the  sinusoidal  character  of  the  signal  waveform, 
will  cause  the  saturation  level  of  the  filter  output  to  be  approximately  equal  to 
that  of  the  filter  input.  Hence,  provided  that  no  right  shifting  is  performed  on 
the  step  3  and  the  step  5  inputs  prior  to  the  filtering,  criterion  (1)  will  be  satis¬ 
fied  for  all  steps  up  to  and  including  step  5,  if  it  is  fulfilled  for  steps  2  and  4  as 
discussed  above.  Nonetheless,  to  provide  a  margin  of  "safety"  for  the  step  3 
filtering,  it  is  prudent  to  take  y  =  1  (i.e.,  to  right -shift  the  step  2  input  by  one 
bit)  even  when  K  =  19.  The  scaling  for  step  6  will  be  discussed  separately. 

Next,  criterion  (2)  of  Section  B.4  must  be  applied.  Particularly,  the  integer 
y  defined  above  must  be  so  chosen  that  the  effect  of  round-off  in  the  final  step  5 
output  is  minimized.  To  check  whether  criterion  (2)  is  satisfied,  the  loss  in 
signal -to -noise  ratio  of  the  step  5  (i.e.,  the  filtered  LASA  beam)  output  must  be 
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calculated  using  equation  (2)  given  in  Section  B.4  and  taking  an  =  0.03  nm.  Thus, 
it  is  necessary  to  calculate  the  bias,  m,  and  the  standard  deviation,  cr  ,  of  the 
round-off  errors  generated  in  the  step  5  output.  This  can  be  done  by  calculating, 
for  each  of  the  steps  2,  3,  4  and  5  in  Figure  B-l,  the  quantum  level,  round-off 
bias  and  round-off  variance  for  the  input  and  the  output.  The  description  of  the 
beamformer  arithmetic  given  in  Section  B.3  will  furnish  the  ground  rules  for 
these  calculations. 

As  an  example,  consider  the  step  3  filtering.  Because  of  the  right-shift  of 
the  step  2  input  by  y  bits,  the  step  3  input  has  the  following  characteristics: 


quantum  level  =  2^Q/K 


round-off  bias  =  a(y)  2y  -  *Q 


round-off  variance  =  a(y) 


Here,  a(y)  =  0,  when  y  =  0  (to  indicate  absence  of  bias  errors  in  the  step  2  output 
when  no  right-shift  is  applied  to  the  step  2  input);  while  a(y)  =  1,  when  y  >  0. 
Therefore,  the  step  3  output  has  the  following  characteristics: 

quantum  level  =  2^Q/K 

round-off  bias  =  ^a(y)A  +  ^2^ 

•8(> 

where:  A,  B,  C,  D  =  the  amplification  factors  of  the  step  3  filters. 

These  factors  are  defined  in  Section  B.2.  The  factor  C  refers  to  the  propagation 
of  random  errors  which  are  statistically  uncorrelated  from  one  sampling  period 
to  another. 


round-off  variance 


a(y)C 
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The  same  technique  can  be  used  for  steps  4  and  5.  Remember  that  a  total 

right-shift  of  5  bits  is  required  prior  to  the  step  4  sum-and-delay  process. 

Hence,  a  right-shift  of  5-y  bits  must  be  applied  to  the  step  4  input.  In  this 

manner,  the  following  expressions  for  the  quantum  level,  the  round-off  bias,  m, 

2 

and  the  round-off  variance,  ,  in  the  step  5  output  are  obtained: 


quantum  level  =  2JQ/KM 


(3) 


m  = 


{°’i  (*«  0  +f)  **6r'I>  +  (cx  f5'1’} 


3KM 


(4) 

(5) 


j|C 

Here,  A^,  Bj,  C^,  C^,  are  the  amplification  factors  of  the  step  5  filter. 

is  the  factor  for  the  propagation  of  that  part  of  the  random  input  error  which  is 

* 

uncorrelated  from  sample  to  sample,  while  applies  to  the  part  of  the  random 

error  which  is  correlated  from  sample  to  sample.  As  explained  in  Section  B.3, 

* 

for  the  filter  considered  here  we  have  C1  <  l.  These  expressions  are  valid  for 

2' 

both  configurations.  Note  that  m  and  crQ  are  minimized  by  choosing  y  as  small 
as  possible.  Hence,  for  minimum  effect  of  round-off  error,  the  amount  of  right- 
shift  of  the  step  2  input  should  also  be  minimum.  Thus,  one  should  select  y  =  0, 
when  K  =  19  and  y  =  1,  when  K  =  25. 

We  shall  now  discuss  each  configurations  separately.  For  both  configurations 
we  take  K  =  19,  M  =  21,  Q  =  0.028 nm,  <rn  =0.03  nm.  For  both  configurations  we 
shall  examine  the  case  where  y  =  0  as  well  as  that  where  y  =  1  (as  pointed  out 
earlier,  y  =  1  is  a  reasonable  choice  even  for  K  =  19). 


B.5.1  First  Configuration 

The  step  3  filter  is  a  sixth  degree  Butterworth  filter  with  the  pass  band  from 
0.8  Hz  to  1.6  Hz.  Its  amplification  factors  are  A  =  0,  B  =  13.7,  C  =  0.17  and 
D  =  560.  Because  the  step  5  filter  is  omitted,  we  have  A^  =  =  1  and 

Bj  =  =  0.  Hence,  by  applying  the  equations  (2),  (3),  (4)  and  (5),  the  results  in 

Table  B-3  are  obtained.  The  quantum  level,  m  and  o^,  are  in  nanometers  (nm). 
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Table  B-3.  Performance  of  First  Configuration 


y 

quantum  level 

m 

a 

e 

loss  (db) 

0 

0.00224 

0.0337 

0.0037 

0.145 

i 

0.00224 

0.0438 

0.0053 

0.422 

Although  the  loss  is  acceptable  for  y  =  1  as  well  as  y  =  0,  note  that  it  is  appreci¬ 
able  for  y  =  1.  The  main  reason  lies  in  the  large  value  of  the  round-off  bias,  m, 
which  exceeds  an  in  both  cases.  Fortunately,  this  bias  is  perfectly  predictable. 
Therefore,  it  can  be  reduced  to  less  than  the  quantum  level  by  adding  to  every 
LASA  beam  an  extra  constant  channel  whose  value  is  the  opposite  of  m.  By  doing 
this,  the  losses  for  y  =  0  and  y  =  1  are  reduced  to  0.065  and  0.145  db,  respectively. 

The  system  performance  is  rather  sensitive  to  the  step  3  filter  choice.  For 
example,  take  y  =  1  and  assume  the  step  3  filter  to  be  a  sixth  degree  Butterworth 
filter  with  a  pass  band  from  0.7  Hz  to  1.4  Hz  (running  at  10  samples  per  second). 

Then  m  =  0.065  and  cr  =  0.009 nm.  The  loss  in  performance  is  2.2  db.  This  loss 

6 

is  excessive,  and  it  is  mandatory  to  compensate  for  the  large  but  predictable 
value  of  m,  using  the  technique  just  explained.  In  this  way,  the  loss  will  be 
reduced  to  0.4  db. 

B.5.2  Second  Configuration 

The  step  3  filter  is  a  low -pass  Butterworth  filter  of  sixth  degree,  with  3-db 
cut-off  at  2.0  Hz  and  running  at  10  samples  per  second.  Its  amplification  factors 
are  A  =  1,  B  =  1.3,  C  =  0.4  and  D  =  1.4.  The  step  5  filter  is  a  sixth  degree  But¬ 
terworth  filter  with  3-db  pass  band  from  0.7  Hz  to  1.4  Hz  and  running  at  five 
samples  per  second.  According  to  Table  B-l,  its  amplification  factors  are 
A.^  =  0,  B^  =  0.8,  Cj  =  0.3  and  D.^  =  5.  The  factor  C^  =  1,  thus  overestimating 
cre  as  well  as  the  loss  in  step  5  output  signal -to-noise  ratio. 

For  y  =  1,  it  is  found  that  m  =  0.0009  and  =  0.0022  nm.  Thus,  the  perfor¬ 
mance  loss  is  0.023  db.  This  loss  is  insignificant  and  will  be  even  less  for  y  =  0. 
Thus,  in  contrast  to  the  first  configuration,  it  is  not  necessary  to  compensate  the 
LASA  beam  outputs  for  round-off  when  using  the  second  configuration. 
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It  is  particularly  instructive  to  compare  the  loss  of  0.023  db  with  the  loss 
of  2.2  db  calculated  above  for  the  same  bandpass  filter,  but  running  it  at  10  sam¬ 
ples  per  second  instead  of  five  samples  per  second,  (see  Section  B.2,  Table  B-l). 

Before  turning  our  attention  to  the  round-off  errors  in  step  6  (see  Figure 
B-l),  it  is  pointed  out  that  a  loss  in  performance  will  occur  as  a  result  of  the 
discrete  phasing  of  the  beamforming  process  in  steps  2  and  4.  This  phasing  is 
implemented  with  a  time  resolution  corresponding  to  a  rate  of  10  samples  per 
second.  The  loss  in  signal  power  is  proportional  to  the  square  of  the  dominant 
signal  frequency  and  is  0.14  db  at  1  Hz  and  0.58  db  at  2  Hz. 

The  results  obtained  so  far  in  this  section  indicate  that,  at  the  price  of  an 
acceptable  loss  in  performance,  the  beamformer  will  produce  an  unsaturated 
waveform  for  signal  amplitudes  up  to  50  nm.  These  waveforms  can  be  used  for 
accurate  event  magnitude  determination,  timing  and  localization  and  for  general 
waveform  analysis. 

To  complete  this  section,  consider  the  rectify -integrate  process  (step  6  in 
Figure  B-l).  In  this  step,  a  short-term  and  a  long-term  average  of  the  rectified 
LASA  beam  output  are  calculated  for  each  LASA  beam  output.  These  two  aver¬ 
ages  are  compared  with  each  other  for  the  purpose  of  event  detection.  We  shall 
not  concern  ourselves  with  the  thresholding  process  but  confine  ourselves  to  the 
calculation  of  the  two  averages. 

We  begin  with  the  short-term  average.  Overflow  can  occur  when  signal  is 
present.  This  will  not  degrade  the  detection  capability,  but  will  affect  the  event 
magnitude,  timing  and  localization  calculations.  However,  these  calculations 
can  then  be  made  using  the  signal  waveform  without  post -detection  averaging. 

As  a  consequence,  overflow  of  the  short-term  rectified  signal  average  does  not 
lower  the  effective  signal  saturation  level  of  the  beamformer. 

No  round-off  error  need  be  incurred  for  the  short-term  average.  For 
example,  straightforward  addition  could  be  used,  always  subtracting  the  oldest 
rectified  sample  while  adding  the  most  recent  rectified  sample.  For  a  rate  of 
five  samples  per  second  and  an  averaging  interval  of  2.5  seconds  this  method 
appears  quite  feasible. 
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This  brings  us  to  the  effect  of  the  post-detection  sampling  rate  of  five  sam¬ 
ples  per  second.  The  noise  standard  deviation  of  the  discrete  sampled  average 
is  always  larger  than  that  of  the  corresponding  analog  average.  However,  the 
loss  in  performance  incurred  is  insignificant  for  the  LASA  beams,  being  on  the 
order  of  only  0.03  db.  In  addition  to  this  loss,  there  is  also  a  difference  between 
the  sampled  and  the  analog  averaged  signal  power.  Sometimes,  the  sampled 
average  is  higher  than  the  analog  average,  and  sometimes,  it  is  lower,  depend¬ 
ing  on  the  location  of  the  signal  onset  relative  to  the  nearest  sampling  point.  The 
term  "signal  power  average"  means  the  maximum  signal  power  average  obtain¬ 
able  for  a  post -detection  integration  interval  of  fixed  duration  but  with  an  arbi¬ 
trary  beginning  point  for  the  analog  case,  and  an  arbitrary  first  sampling  instant 
for  the  sampled  case.  A  quantity  of  interest  is  the  maximum  possible  loss  in 
signal  power  due  to  sampling  (measured  in  db).  This  quantity  not  only  depends 
on  the  post-detection  integration  duration,  T,  and  on  the  sampling  interval,  AT, 
but  also  on  the  signal  waveform.  For  this  reason,  it  can  be  determined  reliably 
only  by  experimental  means.  Nonetheless,  it  is  of  some  interest  to  obtain  a 
rough  estimate.  Assuming  the  signal  to  be  a  sinusoid  of  frequency  f,  it  can  be 
readily  shown  that  the  maximum  loss  is  given  by: 


db 


ForT  =  2  sec  and  AT  =  0.1  sec  (for  a  rate  of  10  samples  per  second),  the  maxi¬ 
mum  loss  is  less  than  0.1  db  when  f  is  less  than  2  Hz.  Thus,  the  sampling  loss 
is  insignificant  for  a  rate  of  10  samples  per  second.  Next,  assume  AT  =  0.2 
sec.  Then,  the  maximum  loss  is  always  less  than  0.2  db  when  f  is  less  than 
1  Hz.  It  is  less  than  0.38  db  when  f  is  less  than  1.5  Hz,  and  less  than  0.78  db 
when  f  is  less  than  2  Hz.  Because  we  are  dealing  here  with  the  "worst  possible" 
situation,  it  would  seem  that  the  losses  are  certainly  acceptable  when  f  is  less 
than  1.5  Hz  (which  is  the  case  for  the  majority  of  the  seismic  events  observed). 
For  a  signal  frequency  of  2  Hz,  the  loss  can  be  as  high  as  0.78  db.  However, 
even  here  it  must  be  emphasized  that  this  is  "worst  case."  For  50  percent  of 
such  signals,  depending  on  the  signal  onset,  the  loss  is  no  more  than  0.2  db. 
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Because  the  computational  load  is  significantly  reduced  when  forming  LASA  beams 
at  five  samples  per  second  rather  than  10  samples  per  second,  it  appears  worth¬ 
while  to  accept  the  performance  losses  associated  with  the  rate  of  five  samples 
per  second  (we  emphasize  that  the  LASA  beam  phasing  must  have  the  resolution 
consistent  with  a  rate  of  10  samples  per  second). 

Next,  consider  the  long-term  average.  It  is  primarily  intended  to  furnish 
the  noise  background  estimate  relative  to  which  the  threshold  operations  must 
be  performed.  The  long-term  average  will  be  calculated  by  feeding  the  rectified 
LASA  beam  output  through  an  exponential  integrating  filter  (i.e.,  y  (m+ 1)  = 
y  (m)  +  A  [x  (m)  -  y  (m)  ] ,  where  y(m)  denotes  the  long-term  average  while  x(m) 
denotes  the  short-term  average  of  the  rectified  LASA  beam  output).  As  noted 
before,  x(m)  can  be  accumulated  without  round-off  error  by  direct  summation. 
Scaling  must  be  implemented  in  such  a  manner  that  neither  overflow  nor  severe 
round-off  errors  occur  in  y(m)  for  integration  times,  t,  on  the  order  of  200  sec¬ 
onds  (t  and  the  filter  coefficient  A  are  related  by  tAf  =  1  where  f  is  the  sampling 

s  s 

rate  for  x(m)). 

To  determine  whether  these  scaling  requirements  can  be  met,  note  that  the 
mean  value  of  y(m)  is  approximately  cr^  while  its  standard  deviation  is  about 
an/T[w:  Here  cr^  designates  the  standard  deviation  of  the  total  (i.e.,  seismic 
plus  round-off)  noise  in  the  step  5  output  while  W  is  its  bandwidth.  A  typical 
value  of  W  is  one  Hz.  From  Section  B.4, 0.03 nm  is  a  reasonable  minimum  value 
of  crn>  Thus,  the  minimum  standard  deviation  for  y(m)  comes  to  0.002  nm,  taking 
t  =  200  seconds. 

The  round-off  error  in  y(m)  has  a  bias  of  Qtf  /2  and  a  standard  deviation  of 

m  s 

Q  1  tf  /24,  where  Q  is  the  quantum  level  of  y(m).  It  is  reasonable  to  demand  a 
s 

scaling  for  which  this  bias  is  small  compared  to  0.03 nm, while  the  round-off 

standard  deviation  is  small  compared  with0.002nm.  Taking  f  =  0.5  samples 

s 

per  second  (corresponding  to  short-term  averaging  over  two  seconds)  and 

Q=10"^nm,  the  round-off  bias  comes  to0.005nmand  the  round-off  standard 

deviation  to  0.00014 nm,  the  scaling  requirements  are  readily  satisfied.  Because 

the  bias  is  predictable,  it  can  be  removed  if  desired. 

15 

The  saturation  level  is  2  Q  =  3.2  nm.  This  level  is  more  than  adequate, 
i.e.,  it  is  unlikely  for  y(m)  to  exceed  that  level.  This  can  be  seen  as  follows.  The 
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seismometer  noise  level  is  on  the  order  of  10  nm  or  less.  Allowing  only  a  20-db 

total  LASA  gain  (instead  of  the  anticipated  30  db,  see  Section  B.3),  one  obtains 

cr  =1  nm.  This  maximum  noise  level  still  is  three  times  below  the  saturation 
n 

level.  Thus,  the  scaling  requirements  for  the  long-term  average  of  the  rectified 
step  5  output  can  be  readily  met. 
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Appendix  C 


PROCESSING  SYSTEM 


C.l  INTRODUCTION 

In  the  First  Quarterly  Technical  Report*  a  functional  system  description 
was  given  and  certain  operating  parameters  were  addressed.  This  Appendix 
extends  some  of  those  considerations,  and  in  addition,  indicates  methods  of 
implementation  of  some  of  the  techniques  suggested  in  Appendices  A  and  B  of 
this  report. 

Section  C.2  discusses  microcoding,  presents  the  algorithms  and  addresses 
the  precision  capability.  The  experimental  display,  a  breadboard  constructed  to 
assist  in  developing  and  evaluating  automatic  'best  beam"  recognition  charac¬ 
teristics  is  discussed  in  Section  C.3. 

In  Section  C.4,  the  overall  system  is  described  in  operational  terms.  The 
functions  of  the  interface  equipments,  maintenance  and  operation  consoles  and 
displays  are  presented.  This  Appendix  concludes  with  Section  C. 5,  an  examina¬ 
tion  of  system  operational  reliability  in  the  detection  mode  and  indicates  the 
improvements  possible  with  the  repairable  duplex  system  approach. 

C.2  MICROPROGRAMMING 

In  the  First  Quarterly  Technical  Report,*  a  system  approach  was  presented 
for  on-line  seismic  detection  using  LASA.  This  approach  involved  the  use  of  five 
signal  processing  algorithms: 


*IBM,  First  Quarterly  Technical  Report,  "LASA  Signal  Procesing  Simulation 
and  Communications  Study,"  Contract  No.  AF  19(628)-5948,  May  27,  1966. 
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a.  Recursive  filter 

b.  Convolution  filter 

c.  Beam  form 

d.  Rectify  and  integrate 

e.  Threshold 

To  achieve  a  maximum  speed  advantage  these  algorithms  can  be  imple¬ 
mented  using  the  technique  of  microprogramming.  This  technique  allows  the 
algorithms  to  be  used  as  if  each  were  a  normal  machine  instruction. 

During  this  quarter,  a  study  and  definition  of  the  algorithm  specifications 
has  been  ur  derway.  One  of  the  primary  premises  for  this  work  was  to  main¬ 
tain  as  much  generality  as  possible  for  the  application  of  each  instruction,  but 
to  also  minimize  the  instruction  execution  time.  The  results  of  the  work  are 
presented  in  the  following  sections. 

C.2.1  Recirsive  and  Convolution  Filters 

These  .two  algorithms  are  combined  because  of  the  similarity  in  their  solu¬ 
tion.  The  two  equations  are: 

a.  Recursive  Filter 


g..(nAt)  =  2 


b.  Convolution  filter 


i  =1.2 


9  •  •  •  9 


J < 


C-2 


For  both  of  the  above  equations,  all  data  and  coefficients  will  have  a  precision 
of  15  data  bits.  The  internal  summations  will  both  be  carried  to  a  precision 
of  31  data  bits,  and  the  outputs  can  be  selected  as  any  15  out  of  the  31  data  bits. 
However,  intermediate  summations  may  extend  upward  to  38  data  bits  with  no 
arithmetic  overflow  indication,  provided  that  the  final  summation  does  not 
exceed  31  data  bits.  Filter  orders  as  large  as  127,  and  i  and  j  subscripts  as 
large  as  32,767  are  allowed. 

The  input  data  for  both  of  these  processes  (f.,  g. j)  is  stored  in  blocks  called 
threaded  lists.  This  means  that  each  block  of  data  representing  one  time 
sample  is  chained  to  the  block  of  data  for  the  previous  time  sample.  By  adjust¬ 
ing  the  chaining  addresses,  processing  rate  changes  can  be  effected  with  a  one¬ 
time  software  change  and  no  change  in  the  algorithms.  This  technique  also 
allows  processing  to  proceed  in  nonconsecutive  sample  periods  using  data  from 
consecutive  sample  periods. 

C.2.2  Beamforming 

The  beamforming  equation  is 


For  this  operation  the  input  data  will  be  15  data  bits  with  a  result  field  of  15 
data  bits.  However,  intermediate  summations  may  extend  upward  to  22  data 
bits  with  no  arithmetic  overflow  indication,  provided  that  the  final  result  does 
not  exceed  15  data  bits.  This  algorithm  can  process  up  to  32,767  beams,  forming 
each  one  from  up  to  32,767  inputs.  The  r-  delay  values  can  have  a  time  resolu¬ 
tion  equivalent  to  that  of  the  input  data,  regardless  of  the  rate  at  which  beam¬ 
forming  is  performed. 

The  beamforming  algorithm  will  use  the  batched  concept  in  which  five  time 
samples  of  a  particular  beam  are  formed  simultaneously.  This  grouping  reduces 
the  computation  time  by  about  43%  because  the  delay  value  must  be  accessed 
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only  once  instead  of  five  times.  The  value  five  was  chosen  because  it  repre¬ 
sented  the  point  at  which  the  percentage  time  gain  started  to  level  off. 

C.2.3  Rectify /integrate  and  Threshold 

These  two  algorithms  are  treated  together  because  their  function  will  be 
combined  into  one  instruction.  The  equations  to  be  solved  are: 


S-l 

Yj(nAt)  =  2  “J  ^  ^  |  x.((n  -  s)At)  | 
S=0 


(1) 


z.(nAt)  =  2n)  £y.(nAt)J  -  2M  J  £z.(  (n  -  l)At)  |+  z.  ((n  -  l)At)  (2) 


l)]+  zj 1 


Zj (nAt)  =  z.((n  -  l)At)  -  y.((n  -  N)At)  +  y.(nAt) 


(3) 


j  1,  2 J. 


In  this  instruction,  all  input  and  output  data  will  be  15  data  bits.  There  is  no 
provision  for  intermediate  summation  overflows. 

The  above  three  equations  are  used  in  solving  for  a  short-term  integration 
value  to  reflect  the  averaged  signal.  Equation  (1)  is  always  solved  followed  by 
either  equations  (2)  or  (3).  This  allows  a  sliding  "window"  look  at  the  data, 
reducing  the  number  of  times  the  instruction  must  be  executed.  In  solving  a 
long-term  integration  for  noise  level,  equations  (2)  or  (3)  are  used.  The  inputs 
to  these  equations  may  be  either  beam  values  or  the  short-term  value  allowing 
alternate  methods  of  computing  the  long-term  average. 
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The  threshold  level  is  set  by  using  the  equation: 


T .  =  u .  *  z.  (nAt) 


(4) 


After  an  event  has  been  detected,  a  different  threshold  can  be  set  by  chang¬ 
ing  the  value  of  u..  The  z.  (nAt)  used  is  the  result  of  the  long-term  inte- 
1  1 

gration  and  this  threshold  is  compared  to  the  short-term  average.  When 
this  threshold  is  exceeded  for  Q  sample  periods  out  of  Q'  consecutive  sample 
periods,  a  detection  is  declared.  The  present  long-term  integration  value  is 
held  until  the  event  no  longer  exceeds  the  threshold,  at  which  time  long-term 
averaging  is  resumed.  The  output  of  the  threshold  section  is  a  list  of  beam 
parameters  for  beams  which  have  exceeded  the  threshold  criterion. 


C.3  EXPERIMENTAL  DISPLAY 

The  investigation  of  methods  of  displaying  event  beams  is  currently  con¬ 
centrated  on  the  pattern  produced  by  the  LASA  beam  intensities.  Beams  are 
formed  and  positioned  according  to  their  relationship  in  the  inverse  velocity 
plane.  Suggested  techniques  to  isolate  the  best  beam  include  determining  the 
center  of  gravity  and/or  moments  of  the  beam  intensity  patterns,  and  pattern 
recognition.  Evaluation  of  these  techniques  would  be  facilitated  by  a  visual 
presentation  of  the  beam  intensity  pattern  as  a  function  of  time.  This  requires 
a  display  with  intensity  modulation  of  a  matrix  of  points  in  the  x-y  plane.  An 
experimental  display  facility  for  this  study  has  been  constructed.  The  display 
utilizes  a  magnetic  tape  for  data  input  and  display  control  and  a  standard  oscillo¬ 
scope  with  Z-axis  input  as  the  output.  A  special  unit  provides  the  interface 
between  the  magnetic  tape  unit  and  the  oscilloscope.  This  interface  unit  includes 
the  necessary  digital-to-analog  converters  and  their  control  logic.  A  block  diagram 
of  the  display  facility  is  shown  in  Figure  C-l. 

The  display  tape  is  machine  generated  with  a  compatible  tape  drive  unit. 

Each  display  tape  runs  for  approximately  6  minutes  which  may  represent  real 
time,  or  be  time  compressed  or  expanded,  as  desired.  Time  scaling  is  per¬ 
formed  when  the  display  tape  is  written. 
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Figure  C-l  Display  Facility  Block  Diagram 


The  display  has  the  capability  of  generating  16  levels  of  intensity  for  1024  beams  in  a 
square  matrix.  Each  beam  appears  as  a  square,  however,  its  size  is  controllable  to  elim¬ 
inate  dark  areas  between  adjacent  beams .  The  patterns  may  be  generated  at  rates  up  to 
6  0  per  second,  thereby,  eliminating  flicker  due  to  the  pattern  fading  between  writings . 

C.4  OPERATIONAL  AND  MAINTENANCE  CONSIDERATIONS 

Within  this  section,  the  processing  equipment  configuration  and  operator  functions 
and  requirements  are  discussed  relative  to  system  design. 

C.4.1  Data  Acquisition 

Raw  LASA  data  is  transmitted  from  each  subarray  equipment  vault  to  the  Data 
Processing  Center  as  digitally  encoded  serial  characters  with  word  parity  checking. 
This  format  is  not  compatible  with  general-purpose  digital  processing  terminals,  and 
thus,  an  Interface  Subsystem  Equipment  (ISE)  requirement  exists  to  preprocess  the 
data.  Maintenance  and  logistics  support  will  be  most  efficient  if  the  same  technology 
is  utilized  throughout  the  center.  Such  a  concept  dictates  compatible  signal  specifi¬ 
cations,  thereby  ensuring  a  minimum  impact  in  the  event  of  subsequent  functional 
changes  or  design  adjustments. 

The  basic  function  of  the  ISE  is  to  buffer  communication  terminals,  multi¬ 
plex  input  channels,  synchronize  the  message  framing,  decode  the  characters 
utilized  in  transmission,  verify  message  fidelity,  and  translate  the  format  to  be 
compatible  with  the  detection  processor  terminal.  A  secondary  function  is  the 
complement  processing  required  for  a  duplex  communication  capability  to  time 
and  control  each  subarray.  The  principal  design  guideline  should  be  programmed 
flexibility  to  ensure  the  initial  interface  with  the  dynamic  subarray  electronics 
configuration,  to  permit  efficient  subsequent  upgrading  of  the  subarray  channels, 
and  to  provide  latitude  for  the  system  evolution  characteristic  of  R&D  efforts. 

Because  of  anticipated  LASA  concept  extensions,  a  growth  potential  is  recom¬ 
mended  to  optimize  the  LASA  program  design  efforts. 

The  most  stable  requirement  item  is  the  physical  capability  to  interface  with 
commercially  available  digital  terminal  equipment  operating  at  one  of  the  standard 
modular  channel  bandwidths.  Bell  half-group  facilities  are  adequate  to  couple 
each  subarray,  and  have  been  utilized  in  the  present  instrumentation. 
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Supporting  simultaneous  communications  through  a  sufficient  number  of  these 
channels  constitutes  a  rudimentary  ISE  baseline  requirement  which  can  be  detailed 
while  retaining  system  flexibility.  Buffering  asynchronous  input  data,  clocking 
output  data,  and  coordinating  terminal  operation  via  the  available  discrete  monitor 
and  control  functions  are  the  principal  interface  design  considerations. 

Bell  303A10  Data  Sets  operate  at  a  19.2  KHz  clock  rate  which  is  slow  when 
compared  with  present  digital  technology.  Early  channel  multiplexing  is  implied 
if  the  ISE  hardware  is  to  be  efficiently  utilized.  Because  this  is  the  first  unifica¬ 
tion  of  data,  redundancy  must  be  provided  to  retain  the  fail-soft  capability  afforded 
by  the  processor  system  architecture.  The  simplest  method  of  achieving  this 
redundancy  with  a  minimum  design  expenditure  is  to  provide  dual  identical  channels 
coupled  to  the  event  and  the  detection  processors,  respectively,  and  to  provide 
instrumentation  to  a  control  and  coordination  channel  enabling  smooth  transitions 
of  data  flow  in  the  event  of  scheduled  maintenance  or  equipment  malfunction. 

Message  formats  are  periodic  with  a  20  Hz  rate,  thereby,  greatly  simplifying 
the  problem  of  maintaining  synchronous  communication.  Two  terminal  modes 
evolve;  the  first  to  initially  capture  frame  coincidence,  and  the  second  to  Verify 
proper  message  format  interpretation  while  operating  with  a  priori  frame  phase 
intelligence.  Residual  effects  of  transients  are  eliminated  by  continuously  refer¬ 
encing  all  communication  timing  to  one  system  standard;  and  because  observed 
propagation  jitter  varies  relatively  slowly,  mode  throwback  would  be  extremely 
rare. 

Parity  error  control  is  presently  employed  at  the  word  level.  Message  bits 
are  encoded  as  characters  with  suboptimal  entropy.  These  must  be  detected  and 
reduced  to  form  a  binary  message  compatible  with  general-purpose  digital  proc¬ 
essors.  Invalid  characters  may  be  correctable  if  a  minimum  Hamming  distance 
error  is  assumed  and  a  most  probable  character  can  be  recognized. 

Because  a  discrepancy  would  have  been  detected,  and  character  substitution 
is  based  on  supposition,  the  word  should  be  flagged  as  marginal  data.  Such  channel 
creditability  qualifiers  may  also  originate  from  alternate  indicators  such  as  word 
parity,  data  set  AGC  level,  or  keys  from  the  maintenance  console. 
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The  exact  format  translations  that  must  transpire  in  the  ISE  are  a  function 
of  the  installation  configuration  and  the  detailed  design  of  the  electronics  located 
at  the  subarray  vault.  Because  these  factors  are  expected  to  vary  as  the  LASA 
is  updated,  a  meritorious  design  feature  would  be  a  programmable  format  trans¬ 
lator.  The  advantages  of  such  an  interface  concept  should  be  comprehensively 
weighed  when  considering  the  slightly  increased  instrumentation  complexity. 

A  complete  history  of  the  LASA  instrumentation  output  and  mode  should  be 
made  available  to  the  recording  equipment  so  that  all  pertinent  data  is  retained, 
thereby  enabling  the  re-creation  of  array  observations  with  a  fidelity  limited  only 
by  the  field  equipment.  However,  a  variable  magnitude  scaling  and  sample  rate 
for  presentation  to  the  detection  processor  is  desirable  to  permit  the  simulation 
of  degraded  field  instrumentation  so  that  design  trade-off  data  could  be  conveni¬ 
ently  obtained  and/or  verified. 

C.4.2  System  Monitor  Requirements 

Supervision  of  a  system  as  intricate  in  equipment  and  concept  as  the  LASA 
data  distribution  and  processing  complex  must  be  based  upon  timely  data  with  an 
adequate  display  of  the  system  status  and  performance.  Consideration  of  these 
essential  system  monitoring  requirements  at  an  early  stage  of  design  can  result 
in  an  efficient  and  orderly  integration  of  the  equipment  required  for  the  system 
supervisory  function.  Figure  C-2  portrays  a  data  monitor  and  crosscheck  capa¬ 
bility  for  event  detection  which  is  consistent  with  the  envisioned  subarray  elec¬ 
tronics  configuration  and  utilizes  the  inherent  flexibility  of  computer-driven  dis¬ 
play  technology. 

Diagnostics  can  be  performed  on  the  analog  instrumentation  by  inserting  sig¬ 
nals  of  known  characteristics  at  various  points  throughout  the  system.  Fault 
detection  and  localization  can  be  accomplished  by  this  technique.  Channel  cali¬ 
bration  is  facilitated  by  a  Proof-Mass  Torqueing  capability  which  superimposes 
disturbance  forces  at  the  transducer.  When  a  broadband  disturbance  of  known 
waveform  is  applied,  the  amplitude  and  phase  characteristics  of  the  channel  can 
be  estimated  as  a  function  of  frequency.  This  channel  signature  will  be  very 
valuable  in  detecting  component  drift  as  a  part  of  system  diagnostics. 
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Figure  C-2.  Data  Monitor  and  Cross-Check  Capability 
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Localization  of  faults  in  the  subarray  vault  equipment  is  presently  based  on 
signal  substitution  or  superposition  through  mode  controls,  and  on  the  monitoring 
of  diagnostic  variables  such  as  temperature  and  key  voltages.  Because  a  good 
portion  of  this  equipment  is  digital,  the  inclusion  of  conventional  digital  diagnostic 
operations  appears  advisable. 

The  telephone  gear  is  well  instrumented  for  communication  diagnostics,  and 
the  details  can  be  found  in  Bell  literature.  In  essence,  it  is  based  on  the  capability 
to  loop  the  channel  at  numerous  points  and  to  insert  artificial  line  losses  into  the 
channel.  Coordination  of  these  tests  with  the  LASA  system  is  beneficial  to  both 
operations  and  diagnostics. 

After  multiplexing,  a  convenient  monitor  tool  is  the  test  channel  which  serves 
as  a  performance  standard.  Inputs  to  this  channel  may  be  deterministic  or  live 
data.  The  output  derived  from  a  deterministic  input  is  easily  checked,  and  pro¬ 
vides  a  means  to  verify  the  process.  Obviously,  an  on-line  channel  cannot  be 
checked  by  this  technique,  and  a  multiplexed  parallel  channel  utilizing  live  data 
may  be  necessary  if  undisturbed  operation  is  essential  during  diagnostic 
operations. 

General-purpose  digital  computers  are  provided  with  a  comprehensive  diag¬ 
nostic  and  error  control  system.  Internal  parity,  format  checks,  processing 
overflow  detection,  and  the  normal  periodic  off-line  diagnostic  operations  ade¬ 
quately  support  the  maintenance  function.  However,  an  on-line  monitor  capa¬ 
bility  is  desirable  to  present  that  data  necessary  for  system  supervision. 

Data  for  the  LASA  processing  system  is  acquired  from  over  500  seismic 
instruments.  Though  each  instrument  is  not  an  indespensible  source  of  intelli¬ 
gence,  each  should  be  monitored  to  assure  timely  detection  at  a  failure  in  the 
signal  integrity.  A  seismic  activity  monitor  as  a  continuous  input  monitor  is 
suggested  to  create  confidence  in  the  input  data  and  to  illuminate  questionable 
channels  so  that  they  may  be  referred  to  the  maintenance  queue  for  detailed 
performance  analysis.  It  is  inconceivable  that  this  broad  monitor  function  could 
efficiently  display  all  maintenance  and  documentation  provisions,  but  rather  that 
it  should  be  a  comprehensive  indicator  of  the  channel  status. 
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A  requirement  exists  to  simultaneously  display  selected  waveforms  to  por¬ 
tray  the  signal  signature,  to  verify  multi-arrivals  from  a  single  event,  and  to 
serve  as  an  editing  table  for  that  data  which  is  selected  for  hard  copy.  An 
observer  will  be  able  to  verify  proper  system  operation  and  participate  in  the 
refinement  of  arrival  time  interpretation.  Input  parameters ,  amplitude  scaling, 
and  time  scaling  must  be  selectable  by  the  operator.  An  adjustable  time  base 
for  each  trace  is  desirable  to  enable  independent  waveform  sliding  for  efficiency 
analysis  of  multi-arrival  events.  This  implies  a  memory  requirement  sufficient 
to  store  each  significant  return  from  a  single  event.  Cursor  measurements  can 
be  compensated  to  ensure  true  time  readings. 

Much  of  the  processed  data  has  geometric  significance  in  "K-space."  Por¬ 
trayal  of  this  intelligence  requires  a  third  display  dimension  which  can  be  effici¬ 
ently  encoded  as  intensity  modulation  of  a  two  dimension  field.  Because  a  large 
dynamic  signal  level  is  encountered,  some  method  of  data  compression  may  be 
desirable.  Decision  theory  suggests  logarithmic  data  compression,  and  this  tech¬ 
nique  has  found  wide  acceptance  in  radar  display  technology.  A  calibrated  pedestal 
and  a  selectable  display  range  are  two  other  features  which  may  be  interpolated 
from  the  extensive  information  retrieval  experience  accumulated  by  the  radar 
industry. 

A  resolution  of  0.5  db  is  desired  to  conveniently  verify  LASA  beam  patterns. 
Figure  C-3  illustrates  the  logarithm  quantum  as  a  variable  number  of  significant 
binary  digits  are  interpreted.  The  required  display  resolution  does  not  materialize 
when  five  or  fewer  digits  are  considered.  Note  that  a  6-digit  algorithm  would  result 
in  a  2nd  quantum  which  is  nearly  2.5  times  the  6th  quantum,  and  would  produce  an 
irregular  display.  Consistency  within  10  percent  of  the  normal  value  is  obtained 
when  a  7-digit  algorithm  is  employed.  This  design  selection  is  considered  opti¬ 
mum  for  the  LASA  System. 

System  maintenance  and  operation  will  involve  console  features  beyond  these 
monitor  displays.  A  digital  input  and  output  will  be  required  to  insert  and  display 
operational  parameters.  Numerous  requirements  for  discrete  indicators  of  status 
and  discrete  control  of  modes  will  evolve  as  the  study  progresses.  Documentation 
facilities  must  be  carefully  configured  to  efficiently  generate  the  essential  hard-copy 
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Figure  C-3.  Logarithm  Quantum  vs  Significant 
Binary  Digits 
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histories  without  burdening  the  facility  with  excessive  logistic  or  storage  require¬ 
ments.  Much  additional  effort  must  be  applied  in  this  area  to  define  a  smoothly 
functioning  and  well  integrated  facility. 

C.5  RELIABILITY  ANALYSIS 

In  the  First  Quarterly  Technical  Report*  a  system  was  presented.  The  pri¬ 
mary  detection  mode  of  operation  required  the  detection  processor  to  continu¬ 
ally  perform  the  detection  process  and  record  seismometer  tape.  The  event 
processor  would  perform  this  function  when  the  detection  processor  was  inop¬ 
erable. 

Of  interest  is  the  reliability  achieved  in  such  a  configuration.  Consider  the 
system  shown  in  Figure  C-4.  The  system  receives  data  from  21  LASA  subarrays 
for  processing.  For  purposes  of  this  analysis,  a  failure  is  the  loss  of  detection 
processing  capability.  The  minimum  requirement  for  detection  processing  is  a 
continuous  path  through  the  reliability  model  shown  in  Figure  C-5  to  a  minimum 
of  three  tape  drive  units.  This  allows  an  event  to  be  detected  and  an  input  tape 
for  event  processing  to  be  written.  The  reliability  analysis  model  can  be  con¬ 
sidered  as  a  cascade  of  four  independent  submodels,  each  of  which  is  required 
for  detection  processing.  Therefore,  the  reliability  of  the  total  system  model  is 
the  product  of  the  reliabilities  of  the  four  submodels. 

Systems  with  redundant  units  which  can  be  repaired  without  disrupting  the 
system  performance  must  be  analyzed  by  techniques  which  consider  the  effect  of 
repair  on  the  effective  system  functional  reliability.  One  method  of  approaching 
the  problem  is  to  employ  a  Markov  process  describing  the  complex  system  rela¬ 
tions  by  the  transitions  or  relationships  between  discrete  states  defined  as  sys¬ 
tem  parameter  analogs.  In  reliability  analysis,  a  state  of  the  system  would  be 
described  by  the  operating  condition  of  the  subunits.  As  an  example,  a  system 
with  two  units  could  be  described  by  the  four  states  in  Table  C-l. 


*IBM,  First  Quarterly  Technical  Report,  "LASA  Signal  Processing  and 
Communications  Study,"  Contract  No.  AF  19(628)-5948,  May  27,  1966. 
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Figure  C-4.  Initial  LASA  Data  Processing  System 
Signal  Flow  Diagram 
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Figure  C-5.  Reliability  Analysis  Model  for  Initial 
LAS  A  Data  Processing  System 
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Table  C-l.  Two- Unit  Operating  States 


State 

Unit  1 

Unit  2 

Operational 

Inoperative 

Operational 

Inoperative 

1 

X 

X 

2 

X 

X 

3 

X 

X 

4 

X 

X 

Whether  the  system  is  operational  for  each  state  would  depend  upon  the  con¬ 
figuration  of  the  units  (cascaded  or  parallel).  Extending  this  to  a  system  with  n 

th 

number  of  units,  the  total  number  of  states  would  be  two  raised  to  the  n  power. 
However,  only  the  states  defining  an  operational  system  need  be  considered.  All 
of  the  other  states  describing  a  failed  system  can  be  grouped  together  as  one 
inoperative  state.  It  is  assumed  that  two  failures  will  not  occur  at  the  same  instant 
of  time  thus  limiting  transitions  to  branches  between  states  differing  by  a  single 
failure.  Figure  C-6a  shows  the  flow  graph  for  the  four  states  of  a  two-unit  system. 
All  four  states  would  be  of  interest  only  for  the  case  of  two  units  being  in  parallel. 
At  analysis  time  equal  to  zero,  both  units  are  assumed  operational.  The  initial 
value  at  the  node  corresponding  to  state  1  would  be  unity  and  the  initial  values 
at  all  other  nodes  would  be  zero.  The  signal  flow  graph  can  be  manipulated  to 
produce  the  configuration  shown  in  Figure  C-6b  when  both  units  have  the  same 
values  of  MTBF  (Mean  Time  Between  Failure)  and  MTR  (Mean  Time  of  Repair). 

The  potential  at  each  node  represents  the  probability  of  the  system  being  in 
the  state  represented  by  the  node.  Therefore,  one  minus  the  potential  as  a  func¬ 
tion  of  time  at  the  node  representing  the  failed  state  gives  the  probability  of 
successful  system  operation  as  a  function  of  operating  time. 
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a.  Two-Unit  System 


Note:  Figure  C-6a  reduces  to  Figure  C-6b  for  two  units 
with  equal  values  of  MTBF  and  MTR  operating  in 
parallel. 


Figure  C-6.  Reliability  Flow  Graph  for  Two-Unit  System 
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Using  typical  representative  assumed  values  of  MTBF  and  MTR,  the  prob¬ 
ability  of  successful  uninterrupted  processing  as  a  function  of  operating  time 
was  calculated  for  the  system  configuration  shown  in  Figure  C-4.  The  prob¬ 
ability  of  successful  system  operation  for  100  hours  with  and  without  repair 
is  0.958  and  0.830,  respectively,  and  for  200  hours,  0.914  and  0.557,  respec¬ 
tively.  The  marked  improvement  in  system  operation  with  repair  indicates  the 
need  for  both  scheduled  preventive  maintenance  as  well  as  "on-call"  maintenance. 
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Appendix  D 


COMMUNICATION  STUDY 

Data  transmission  represents  a  significant  fraction  of  the  LASA  system 
functions.  The  19.2Kbs  data  rate  from  each  of  the  21  subarrays  generates 
a  serial  403.2  Kbs  data  rate  when  multiplexed  onto  a  single  line.  Due  to 
the  redundant  bit  coding  used  in  the  subarrays,  the  composite  information  rate 
is  one-half  this  data  rate,  or  201.6 Kbs.  In  this  appendix,  the  transmission  of 
LASA  data  from  a  collection  point  to  a  remote  signal  processing  center,  and 
certain  trade-offs  that  may  be  relevant  in  reducing  the  data  transmission  load, 
are  addressed.  In  addition,  some  comments  on  message  framing  and  error 
encoding  techniques  pertinent  to  the  long-line  transmission  techniques  are 
included. 

D.  1  DATA  TRANSMISSION 

To  provide  an  estimate  of  the  costs  involved  in  data  transmission,  a  specific 
case  of  transmitting  LASA  data  from  the  present  Billings,  Montana  facility  to 
a  signal  processing  center  located  in  the  Washington,  D.  C.  area  was  examined. 
Although  several  locations  were  considered,  only  one  is  presented  here  to  exhibit 
method  of  analysis  and  trade-off  variables.  Two  approaches  were  considered. 

The  first  assumed  that  data  must  be  transmitted  automatically  and  in  real  time. 
The  second  assumed  that  the  fixed  delay  inherent  in  shipping  the  data  via  air 
freight  was  acceptable.  In  both  cases, the  costs  involved  in  acquiring  the  data 
have  been  omitted  because  the  added  cost  of  the  tape  writing  operation  is  approxi¬ 
mated  by  the  cost  for  an  added  unit  between  the  modem  and  the  remote  signal 
processing  system. 


D-l 


The  Commercial  Communication  Rates  for  equipment  covering  the  range  of 
bandwidths  of  interest  are  presented  in  Table  D-l  and  the  costs  for  manual  data 
transmission  are  presented  in  Table  D-2.  Using  Tables  D-l  and  D-2,  the  auto¬ 
matic  data  transmission  cost  compares  favorably  to  the  manual  method  for  data 
rates  accommodated  by  both  the  TELPAK  A  and  TELPAK  B.  Without  attaching 
any  economic  significance  to  the  fidelity,  reliability  and  convenience  of  the  auto¬ 
matic  method  over  the  manual  method,  an  additional  trade-off  variable  of  dis¬ 
tance  between  the  remote  processing  site  and  the  central  collection  point  can  be 
examined.  Because  most  of  the  cost  is  a  direct  function  of  distance,  the  cost 
break  point  in  favor  of  TELPAK  C  and  D  exists  for  tariff  mileages  less  than 
1333  and  725,  respectively. 

The  data  rate  at  Billings  would  require  a  TELPAK  D  system,  however,  the 
useful  information  rate  can  be  accommodated  by  a  TELPAK  C.  Because  TELPAK 
B  exists  physically  as  two  TELPAK  A's  operated  in  parallel,  its  utilization  intro¬ 
duces  complicating  control  problems.  Therefore,  efforts  to  reduce  data  rates 
below  TELPAK  C  capability  have  been  aimed  at  fitting  within  a  TELPAK  A 
channel. 

Three  variables  have  been  considered  in  trade  offs  aimed  at  reducing  the 
data  rate : 

a.  Sampling  Rate— The  current  20  Hz  sampling  rate  appears  excessive 
with  respect  to  the  anticipated  signal  frequency  and  further  analysis 
may  verify  that  a  10  Hz  sampling  rate  is  adequate.  This,  combined 
with  the  reduction  in  redundant  bits  permitted  by  polynomial  coding 

for  message  framing  and  error  checking,  allows  use  of  100 Kbs  channels. 

b.  Number  of  Words  per  Subarray— The  32  words  received  from  each 
subarray,  each  sample  period,  presently  contain  only  26  words  of 
interest  in  beamforming.  In  addition,  current  simulation  processing 
casts  doubt  upon  the  usefulness  of  the  inner  ring  of  six  seismometers 
per  subarray,  offering  further  potential  reduction  in  the  number  of 
words  to  be  transmitted. 
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Table  D-l.  Commercial  Communication  Rates 


Kc 

Kbs 

$/mile 

$ /Month 

Two  Terminals 

One-Time 

Installation/2 

Voice  Line 

4 

2.4 

$  2.02  1st  150 
1.717  to  500 
1.616  over  500 

$50 

$20 

TELPAK  A 

48 

50 

15.00 

680 

500 

TELPAK  B 

96 

100 

20.00 

1360 

1000 

TELPAK  C 

240 

230.4 

25.00 

1890 

1000 

TELPAK  D 

480 

500 

45.00 

2600 

1800 

Monthly  Costs,  Billings  to  Washington. 
Tariff  Mileage  =  1631  miles 


Voice  Line 


TELPAK  A 


TELPAK  B 


TELPAK  C 


TELPAK  D 


Mileage 
Terminals 
Total  Monthly 


$2,950 

50 

$3,000 


$24,500 

680 

$25,180 


$32,600 

1,360 

$33,960 


$40,800 

1,890 

$42,690 


$73,395 

2,600 

$75,995 


Note:  A  recent  U.  S.  Appeals  Court  decision  has  ordered  AT&T  Corp.  to  make 
TELPAK  A  and  B  rates  correspond  to  private  line  rates  and  to  increase  TELPAK 
C  and  D  rates.  Major  TELPAK  users  have  asked  for  a  stay  of  this  order  until 
they  can  carry  an  appeal  to  the  Supreme  Court.  Ultimate  effect  of  this  litigation 
on  TELPAK  rates  is  uncertain  at  this  time. 
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Table  D-2.  Data  Transmission  Costs  by  Tape 

(Billings,  Montana  to  Washington,  D.  C.) 


Assumptions 

a.  Operation  is  24  hours  per  day  and  7  days  per  week. 

b.  800  BPI,  3  characters  per  word  format  yielding  160  tapes  per  day. 

c.  All  operations  done  in  parallel  with  present  recording  activity. 

d.  CPU  time  is  available  to  write  second  tape. 

t 

Purchase 

Monthly 

1.  Tape  Control  Unit 

$18,900.1 

2.  Two  Tape  Drives 

60,800  .X 

3 .  Montana— two  operators  per  shift 

$11,700. 

@  $1300  per  man  month 

4.  Air  Freight  @  $53  per  100  lbs. 

10,812. 3 

5.  Washington— two  operators  per  shift 

11,700. 

@  $1300  per  man  month 

6.  Tape  inventory,  1600  reels 

53,600.2 

@  $33 .50  per  reel 

7.  Estimate  Maintenance  Cost 

1,000. 

$133,300. 

$35,212. 

1.  Equipment  catalog  costs . 

2.  Costs  from  IBM  Authorized  Federal  Supply  Schedule  Price  List  dated 
January  17,  1966  for  800  BPI  tested  tape  IBM  P.  N.  339269. 

3.  Air  freight  costs  based  on  4.25  pounds  per  tape  reel  and  $51  round  trip 
air  freight  with  $1  ground  handling  at  each  terminal  for  each  100  pounds 
shipped. 
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c.  Word  Size— Work  is  now  in  progress  to  determine  the  necessity  or 

desirability  of  carrying  the  full  13  bit  plus  sign  quantization  of  the 
data  words.  In  addition,  it  has  been  shown  that,  for  low  signal-to- 
noise  ratio,  beam  forming  with  infinitely  clipped  signals  results  in 
negligible  loss.* 

A  subarray  message  configuration  which  may  be  worthy  of  consideration  uses 
19  words  for  short-period  instruments,  three  words  for  long-period  instruments, 
three  for  environment  sensors,  one  control  word,  and  two  spares.  Using  Figure 
D-l,  it  can  be  shown  that  employing  a  28-word  message  allows  8  bits  per  word 
at  10  Hz  in  a  TELPAK  A. 

D.2  OTHER  CONSIDERATIONS 

For  long-distance  transmission,  it  appears  necessary  to  not  only  provide 
message  framing  but  also  error  detecting.  The  Start  of  Message  word  for  mes¬ 
sage  framing  and  simple  word  parity  of  the  present  system  utilizes  approximately 
10%  of  the  available  data  bits  and  is  not  particularly  effective  for  the  burst  and 
dropout  errors  that  are  characteristic  of  long-line  communication  links .  An 
alternate  approach  would  be  to  utilize  a  polynomial  coding**  as  an  efficient  means 
of  providing  simultaneous  error  coding  and  message  framing.  Assuming  that 
data  is  transmitted  in  1000  bit  messages,  a  28-bit  polynomial  code  has  been 
selected  for  message  framing  and  error  checking.  This  represents  only  2.8% 
redundant  bits  for  both  functions,  providing  a  strong  message  framing  capability 
and  virtually  eliminating  undetected  transmission  errors . 


*V.  C.  Anderson,  "Digital  Array  Phasing,"  Journal  of  the  Acoustical  Society  of 
America,  Vol.  32,  No.  7,  1960. 

*P.  Rudnick,  "Small  Signals  in  the  DIMUS  Array,"  Journal  of  the  Acoustical 
Society  of  America,  Vol.  32,  No.  7,  1960. 

**A.  H.  Frey,  Jr.,  "Message  Framing  and  Error  Control,"  IEEE  Transactions  an 
Military  Electronics,  Vol.  MIL9,  No.  2,  April  1965. 
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Assumptions: 

a.  Available  channel  is  50  kbs,  TELPAK  A. 

b.  Messages  are  1/000  bits  long. 

c.  A  28-bit  polynomial  is  used  for  message  framing 

and  error  detection. 

d.  Sampling  Rate  =  10  Hz. 

Figure  D-l.  Word  Size  vs  Number  of  Words 
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The  specific  28-bit  polynomial  was  generated  in  the  following  manner :  The 
first  two  irreducible  polynomials  of  degree  10  in  Peterson's  Tables  were  multi¬ 
plied  together  to  obtain  a  polynomial  of  degree  20.  This  polynomial  is  then 
3  5  2 

multiplied  by  X  +1  and  X  +  X  +  1  to  obtain  the  resulting  28th  degree  polynomial: 
X28  +  X23  +  X20  +  X19  +  X15  +  X14  +  X11  +  X10  +  X8  +  X4  +  X+l. 


The  polynomial  code  thus  formed  has  a  minimum  distance  of  6,  ensuring 
detection  of  all  error  patterns  comprising  five  or  fewer  erroneous  bits,  regard- 

3 

less  of  location  in  the  message.  The  incorporation  of  the  factor  X+l  (in  X+l) 
ensures  that  all  patterns  with  an  odd  number  of  errors  will  be  detected.  The 
error  detection  capability  can  be  summarized  as  follows: 

a.  All  error  patterns  containing  five  or  fewer  errors  in  the  message  are 
certain  to  be  detected. 

b.  All  single  bursts  of  error  of  duration  less  than  29  bits  are  certain 
to  be  detected. 

c .  Of  the  error  patterns  with  six  or  more  errors ,  all  patterns  with  an  odd 
number  of  errors  will  be  detected. 

d.  Error  patterns  with  six  or  more  bits  in  error  will  be  detected  with 

a  probability  approximately  equal  to  1  minus  the  quantity  2  raised  to 

-28 

the  minus  28th  power  (1-2  ). 

Based  on  the  error  data  given  by  Fontaine  and  Gallager*  for  voice  band 
channels ,  an  undetected  error  for  the  Bell  A-l  modem  would  occur  on  the  average 
of  once  every  18.5  thousand  years.  Although  data  is  not  available  for  the  Bell 
303A  series  of  modems  utilized  with  the  LASA  arrays ,  somewhat  comparable 
results  would  be  anticipated  and  undetected  errors  can  be  expected  to  be  of  no 
consequence. 

Implementation  of  the  encoding  and  decoding  of  such  polynomials  can  be  done 
in  a  general  purpose  computer,  but  not  very  efficiently. 


♦Fontaine  and  Gallager,  "Error  Statistics  and  Coding  for  Binary  Transmission 
Over  Telephone  Circuits"  Proceedings  of  the  IRE,  June  1961. 
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Special  polynomial  encoder-decoders  have  been  built,  however,  that  efficiently 

implement  this  operation.  As  an  example,  the  encoding  and  decoding  logic  for 

2  4  5 

the  polynomial  P(X)=1+X  +X  +X  is  shown  in  Figure  D-2. 

For  encoding,  the  message  enters  the  input  of  the  shift  register  to  generate 
the  five  polynomial  bits  as  indicated  in  the  logic.  The  AND  gate,  A^,  is  closed 
to  allow  the  feedback  through  the  EXCLUSIVE  OR  gates  as  the  five  redundant 
bits  are  generated.  Simultaneously, the  message  is  fed  through  the  five-bit  delay 
to  the  output.  After  all  data  bits  plus  five  zeros  have  been  shifted  in,  the  remainder 
in  the  register  constitutes  the  five  redundant  bits.  AND  gate  A.^  is  then  opened 
and  AND  gate  A ^  closed  and  the  five  redundant  bits  sent  out  at  the  end  of  the  mes¬ 
sage.  Decoding  operates  in  the  same  manner,  with  the  error  check  for  zero  after 
the  received  message  and  the  redundant  bits  have  been  put  through  the  register . 

To  incorporate  message  framing,  a  message  length  buffer  is  added  because 
the  receiving  end  does  not  know  where  the  start  of  a  message  is  until  it  has 
received  the  entire  message  with  its  redundant  bits  and  completed  the  polynomial 
check.  In  effect,  the  register  looks  at  the  last  N  bits  where  N  is  the  total  mes¬ 
sage  length  (data  plus  redundant  bits),  to  determine  if  a  valid,  error-free  message 
exists .  Thus ,  the  data  is  continually  scanned  for  valid  messages  through  a  sliding 
window,  N  bits  long. 

In  this  type  operation,  one  other  addition  must  be  made.  Because  equipment 
is  always  looking  at  the  last  N  bits  to  see  if  it  is  a  valid  message,  it  does  not  want 
to  be  biased  by  earlier  bits .  Therefore ,  after  N  bits  have  been  shifted  through 
the  register,  when  the  N  +  1  bit  is  shifted  in,  the  first  bit  is  fed  back  into  the 
register  to  cancel  its  effect  on  the  running  value  in  the  polynomial  register.  In 
this  way, for  both  encoding  and  decoding,  the  register  always  contains  the  value 
corresponding  to  the  last  N  bits  that  have  gone  through  it.  In  actual  implementa¬ 
tion,  other  features,  such  as  a  nonzero  check  value,  may  be  added  to  guard  against 
the  more  common  catastrophic  failure  conditions. 

The  value  of  the  above  type  of  error  detecting  lies  not  only  in  the  protection 
it  gives  the  signal  processing  system  from  transmission  noise,  but  in  the  strong 
diagnostic  aid  it  provides  for  fault  conditions.  Because  for  all  practical  purposes, 
the  undetected  error  rate  is  zero,  such  a  system  allows  immediate  isolation  of  the 
transmission  line  as  a  fault  source. 
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FF— F I  ip— Flop  Register  Position 


Figure  D-2.  Coding  Logic 
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DATA  ANALYSIS  PROGRAMS 


E.l  INTRODUCTION 

This  appendix  presents  a  general  description  of  the  Data  Analysis  Programs 
being  developed. 

For  the  most  part,  all  the  programs  discussed  are  operational  in  the  sense 
that  a  working  model  of  each  program  is  available.  Most  of  the  programs,  how¬ 
ever,  are  not  complete  and  changes  are  being  made  either  to  increase  their 
capability  in  a  general  sense  or  to  modify  them  to  the  individual  characteristics 
of  a  particular  experiment. 

The  programs  are  generally  written  in  FORTRAN  IV  so  that  they  are  easily 
adaptable  for  processing  on  different  machines.  The  majority  of  the  programs 
are  currently  being  executed  on  an  IBM  7090  although  some  are  specifically 
written  for  the  IBM  System/360.  The  input  routines  for  the  LASA  Tape  Edit  and 
Subarray  Beamformer  programs  are  written  in  IBM  7090  MAP  language  and 
mast  be  rewritten  in  System/360  Assembler  Language  prior  to  conducting 
experiments  on  the  System/360. 

E.2  LASA  TAPE  EDIT 

E.2.1  Purpose 

The  LASA  Tape  Edit  program  converts  LASA  raw  data  tapes  into  FORTRAN 
compatible  binary  tapes.  Any  or  all  651  channels  or  any  portion  thereof  may  be 
selected  for  editing. 
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E.2.2  Description 


E.2.2.1  Input 

Program  input  consists  of  edit  processing  control  cards  and  from  one  to 
three  LASA  raw  data  tapes. 

The  control  cards  are  used  to  indicate  to  the  program  (a)  the  number  of 
records  to  be  skipped  before  initiating  the  editing  process,  (b)  the  number  of 
records  to  be  edited,  (c)  the  number  of  channels  to  be  edited,  (d)  the  number  of 
LASA  raw  data  tapes  to  be  mounted  and  their  respective  names  and  reel  numbers, 
and  (e)  the  channel  number  of  each  channel  to  be  edited. 

Each  LASA  raw  data  tape  is  a  low-density  (556  BPI)  tape  with  each  data 
word  being  18  bits  in  length  with  extended  sign  and  odd  parity.  Negative  numbers 
are  in  two's  complement  notation.  A  record  on  a  raw  data  tape  consists  of :  (a)  a 
four-word  header  containing  the  date-time  data,  (b)  two  frames  each  containing 
651  words  representing  one  data  sample  from  each  of  the  651  channels  per  time 
interval  of  50  milliseconds,  and  (c)  a  seismometer  status  trailer  when  required. 
All  trailer  data  is  ignored  during  the  editing  process. 

E.2.2. 2  Computation 

Each  record  on  the  raw  data  tapes  is  unpacked  and  placed  in  individual 
IBM  7090  words.  The  data  channel  values  (seismometers)  are  converted  to  true 
floating-point  notation.  The  data  is  then  arranged  in  memory  for  subsequent 
recording. 

E.2.2. 3  Output 

Program  output  consists  of  an  off-line  printout  and  as  many  FORTRAN 
compatible  LASA  edit  tapes  as  are  required  during  the  editing  process. 

The  off-line  printout  lists  all  the  information  supplied  on  the  input  control 
cards  as  indicated  in  Section  E.2.2.1.  It  also  contains  the  number  of  logical 
records  on  the  output  edit  tapes,  raw  data  tape  record  numbers  where  read 
checks  occurred  (if  any),  the  locations  of  unexpected  end  of  files  on  the  raw 
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data  tapes  encountered  prior  to  the  completion  of  the  specified  number  of 
records  to  be  edited,  and  data  from  the  header  record  written  on  the  first  LASA 
edit  tape  as  described  below. 

Each  LASA  Edit  tape  is  a  high-density  (800  BPI)  tape  with  each  data  word 
being  36  bits  in  length.  Each  logical  record  on  LASA  edit  tape  will  consist  of 
from  one  to  651  data  words.  Because  these  logical  records  are  written  in 
FORTRAN,  each  will  be  subdivided  into  physical  records  of  no  more  than  256 
words.  The  first  logical  record  is  a  header  record  which  contains  in  consecu¬ 
tive  word  order  the  year,  day,  hour,  minute,  second,  and  millisecond  of  the  first 
data  sample  record  which  follows  the  header,  and  the  number  of  edited  channels 
and  their  respective  channel  numbers.  The  header  words  are  in  fixed-point 
form.  Each  logical  record  which  follows  the  header  record  consists  of  a 
floating-point  data  sample  of  all  the  channels  selected  for  editing.  As  raw  data 
tapes  are  recorded,  four  data  samples  are  lost  during  reel  switching.  To 
compensate  for  these  lost  records,  four  logical  records  of  the  last  available 
data  will  be  recorded  when  an  end  of  file  is  reached  on  all  but  the  last  raw  data 
tape . 


E.2.3  Program  Interaction 

The  edited  tapes  generated  by  this  program  form  the  input  for  the  Subarray 
Beamformer  (SBF)  program. 

El  .2.4  Program  Restrictions 

The  program  expects  only  one  end  of  file  per  LASA  raw  data  tape.  Because 
IASA  raw  data  tapes  do  not  contain  a  year  in  their  header,  an  NYEAR  =  card 
must  be  inserted  in  the  "Read  I  Subroutine ."  The  present  deck  has  a  NYEAR  = 
1966. 

E.2.5  Comments 

The  LASA  Tape  Edit  program  is  designed  for  the  IBM  7090  computer  under 
I3SYS  supervision  and  is  written  in  FORTRAN  IV  and  MAP.  All  reading  is  done 
by  Library  IOCS  and  writing  by  FORTRAN  IV. 
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E.3  SIGN AL-TO -NOISE  RATIO  CONTROL 


E.3.1  Purpose 

This  program  constructs  a  LASA  data  tape  containing  a  specified  number 
of  noise  samples  followed  by  a  specified  number  of  signal  samples,  each  of 
which  has  been  multiplied  by  a  unique  scaling  factor  and  added  to  a  noise  sample. 
The  overall  effect  is  the  generation  of  a  magnetic  tape  containing  an  event  in 
which  the  signal -to -noise  ratio  can  be  controlled.  This  tape  would  then  be  avail¬ 
able  for  processing  through  pertinent  data  analysis  programs. 

E.3.2  Description 

E.3. 2.1  Input 

The  program  reads  a  LASA  tape  in  edited  form  and  extracts  segments  of 
noise  samples  and  signal  samples.  Simple  control  cards  provide  the  capability  to 
precisely  choose  the  portions  of  noise  and  signal  samples  to  be  used  by  the  pro¬ 
gram  in  generating  a  pseudo-LASA  event  tape  with  a  known  signal-to-noise 
ratio.  Scaling  factors  for  each  seismometer  present  on  the  input  tape  are  pro¬ 
vided  by  card  input. 

E. 3.2.2  Computation 

A  pseudo-event  tape  may  consist  of  more  than  one  reel  (4  minutes  of  LASA 
seismometer  data)  and  may  use  more  than  one  reel  of  input  tape  for  its 
generation. 

The  program  has  the  following  characteristics: 

a.  Each  channel  (seismometer)  must  have  a  distinct  scaling  factor.  If 
there  are  more  channels  than  scaling  factors,  the  remaining  channels 
will  be  scaled  by  a  factor  of  1.0. 

b.  A  designated  number  of  time  samples  may  be  skipped  prior  to  taking 
the  noise  sample. 
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c.  A  designated  number  of  records  may  be  skipped  between  the  end  of 
the  noise  sample  and  the  beginning  of  the  signal. 

d.  The  number  of  noise  samples  must  be  greater  than  the  number  of 
signal  samples.  The  difference  will  be  the  number  of  noise  samples 
at  the  beginning  of  the  output  tape. 

E. 3.2.3  Output 

Output  tape(s)  are  written  in  the  same  format  as  the  input  tape(s)  and 
contain  the  same  number  of  channels.  The  first  X  number  of  records  are 
noise  samples,  and  the  rest  of  the  records  are  scaled  signal  plus  noise. 

E.3.3  Program  Interaction 

The  input  and  output  tapes  are  in  the  LASA  Tape  Edit  format.  Therefore , 
signal -to-noise  ratio  can  be  adjusted  on  any  tapes  of  that  format. 

E.3.4  Program  Restrictions 

The  only  restriction  is  that  enough  data  must  be  present  on  the  input  tape 
to  accomplish  the  production  of  the  signal-to-noise  ratio  output. 

E.3.5  Comments 

The  program  is  written  in  FORTRAN  IV  and  MAP  for  the  IBM  7090. 

E.4  INVERSE  VELOCITY  SPACE  MAPPING 
E.4.1  Purpose 

This  program  prepares  inverse  velocity  space  maps  of  the  world  land  areas 
and  seismic  zones  as  well  as  lines  of  constant  geocentric  latitude  and  longitude. 
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E  .4.2  Description 


E.4.2.1  Input 

Two  forms  of  input  are  acceptable:  (a)  pairs  of  latitude  and  longitude  to  be 
mapped  into  inverse  velocity  coordinates,  or  (b)  initial  latitude  and  longitude 
coordinates  and  an  incremental  value  for  each  to  generate  coordinate  lines  for 
inverse  velocity  space  mapping. 

E.4.2.2  Computation 

The  transformations  consist  of  calculating  range  and  azimuth  by  means  of 
spherical  trigonometry  and  then  converting  range  to  horizontal  phase  velocity. 
For  more  detailed  discussion  of  the  computation  involved  in  this  program,  see 
Appendix  A . 

E.4.2.3  Output 

The  output  of  this  program  consists  of  a  printout,  for  each  coordinate,  of 
range,  the  inverse  phase  velocity, and  azimuth;  and  a  plot  produced  for  the 
CALCOMP  plotter. 

E.4.3  Program  Interaction 
None. 

E.4.4  Program  Restrictions 
None. 

E.4.5  Comments 

The  program  is  written  for  the  IBM  7090  in  FORTRAN  IV  using  CALCOMP 
plotter  subroutines. 
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E.5  FILTER  PROGRAMS 


E.5.1  Coefficient  Generator  Program 
E.5. 1.1  Purpose 

This  program  generates  recursive  filter  coefficients  for  a  given  low-pass, 
high -pass,  band-pass  or  band-reject  Butterworth  or  Chebychev  filter  with  the 
option  of  calculating  a  frequency  response,  an  impulse  response  of  both  the  filter 
transfer  function  and  its  denominator,  and  corresponding  bias  and  variance 
coefficients. 

E.5. 1.2  Description 

E.5. 1.2.1  Input 

The  desired  filter  variables  (type,  kind,  order)  must  be  specified  along 
with  sampling  rate  and  center  and/or  cut-off  frequencies. 

E.5.1.2.2  Output 

This  program  gives  as  printout  the  a  and  b  filter  coefficients  along  with  a 
frequency  response. 

E.5. 1.3  Program  Interaction 

Filter  coefficients  generated  by  this  program  will  be  used  in  the  recursive 
filter  in  the  Subarray  Beamformer  (SBF)  program. 

E.5. 1.4  Program  Restrictions 

Subject  to  precision  requirements,  the  program  is  limited  to  order  20  for 
low-pass  and  high-pass  filters  and  limited  to  order  10  for  band-pass  and  band- 
reject  filters. 
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E.5.1.5  Comments 


This  program  is  written  in  FORTRAN  IV  for  7090.  The  equations  used  in  the 
computations  are  not  finalized  and  are  not  presented  herein. 

E.5.2  Frequency  Response  Program 

E.  5.2.1  Purpose 

This  program  generates  a  frequency  response  for  a  set  of  input  filter 
coefficients. 


E.5.2 .2  Description 


E.5.2. 2.1  Input 

The  a  and  b  filter  coefficients  are  required  as  input. 


E. 5.2.2. 2  Computation 

The  magnitude  of  the  loss  in  db  and  the  phase  angle  are  calculated  for  a 

range  of  frequencies, 

k  k 


i  +  y  a 
0  Li  n 


(  )  a  sii 

L  n 


H  = 


(a^  +  )  a_  cos  nwT)  -  j  (  >  an  sin  nwT) 

~n=l _ n=l _ 

k  k 

bn  cos  nwT)  -  j  bR  sin  nwT) 


n=l 


n=l 


H 


P  -3Q 

R  -  jS 


where: 


f. 

wT  = 

xs 


a  and  b 
n  n 


the  filter  coefficients 

the  order  of  the  filter 

the  frequency  being  investigated 

the  sampling  rate . 
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The  magnitude  of  the  loss  in  db  is  expressed  as 

l»l  -  101°*iofr^r 

and  the  phase  angle  is  expressed  as 


<  H 


arctan 


PS  -QR 

PR  +  QS 


E. 5.2. 2. 3  Output 

The  printout  consists  of  a  range  of  frequencies  with  the  corresponding  db 
loss  and  phase  angle  for  each  of  these  frequencies. 

E.5.2.3  Program  Interaction 
None. 


E.5.2.4  Program  Restrictions 

The  program  has  been  restricted  to  a  maximum  of  20  filter  coefficients. 
E.5.2.5  Comments 

The  program  was  written  in  FORTRAN  IV  for  the  IBM  7090. 

E.6  SUBARRAY  BEAMFORMER  (SABF) 

E.6.1  Purpose 

This  program  produces  a  variable  number  of  subarray  beams  for  each 
subarray  considered.  In  addition,  the  option  to  pass  the  beams  through  a  digital 
filter  is  also  available.  The  purpose  of  the  program  is  to  produce  unfiltered  and 
filtered  subarray  beams  which  will  be  used  for  further  processing. 
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E.6.2  Description 


E.6.2.1  Input 

LASA  edited  tape(s)  is  the  main  input  to  the  SABF  program.  This  input 
represents  the  data  from  a  selected  number  of  seismometers  for  a  selected 
period  of  time.  The  number  of  seismometer  channels  recorded  on  the  tape 
dictates  the  number  of  elements  in  each  beam  sum.  An  equal  number  of 
seismometers  from  each  subarray  must,  however,  be  present  on  the  tape. 

In  addition  to  the  seismometer  data,  digital  filter  parameters  and  time 
delays  form  the  other  input  data. 


E.6.2. 2  Computation 

There  are  three  main  computations  present  in  SABF:  beamforming, 
filtering,  and  thresholding. 


E. 6.2.2. 1  Beamforming 

The  following  equation  represents  the  computation  involved  in  the  process 
of  beamforming  in  SABF: 

I 

vt)=fl  v'-v 

i=l 


where: 

I 


3 

k  = 
t 


f..  = 

ij 

v(t)  = 


number  of  seismometers 
•th  , 

j  subarray 
k**1  beam 

current  time  sample 

time  delay  for  the  i^  seismometer  in  the  j*  subarray,  for  the 
k**1  beam 

LASA  raw  data  from  the  i^seismometer  in  the  j^1  subarray 

th  th 

resultant  subarray  beam  for  the  k  beam  in  the  j  n  subarray. 
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As  indicated  in  the  equation,  each  beam  formed  represents  a  delayed  sum  of  a 
set  of  seismometers  from  a  particular  subarray.  This  process  is  repeated, 
using  a  different  set  of  time  delays,  for  each  beam  formed  in  that  subarray, 
and  then  repeated  for  the  total  number  of  subarrays  present.  The  number  of 
beams  formed,  then,  is  j  times  k,  with  each  beam  formed  by  summing  I  properly 
delayed  seismometers.  The  variables  I,  j,  k,  as  well  as  the  data,  f„,  and  time 
delays,  n.^,  are  provided  as  input  to  this  program. 


E. 6. 2.2 .2  Recursive  Filtering 

The  following  equation  represents  the  computation  involved  in  the  process 
of  recursive  filtering  in  SABF: 


P  P 

Gjk  <l>  ■  Y  ajk  ^  '  Bjk (t_p)  •  Y  V  <e>  •  Gjk  <‘-P> 

p=o  p=l 


where: 

P 


j 

k 

t 

ajkte) 

bjk(P) 

Bjk(t-p) 

Gjk(t-p) 


the  order  of  the  filter 
•th  , 

j  subarray 
k**1  beam 

current  time  sample 
filter  coefficients 
filter  coefficients 

th  th 

the  km  beam  formed  in  the  j  n  subarray  at  sample  time  (t-p) 

th  th 

the  k  filtered  beam  formed  in  the  j  subarray  at 

sample  time  (t-p) 


As  the  equation  indicates,  the  filtering  is  performed  on  the  beams  formed  in 
SABF.  Each  beam  formed  is  passed  through  the  same  digital  filter  and  a  resul¬ 
tant  filtered  beam  value,  Gjk(t)>  for  the  current  time  period  is  produced.  The 
variables  P,  j,  k,  a,  and  b  are  provided  as  input  to  this  program. 
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E. 6.2.2. 3  Thresholding 


Each  beam  and  filtered  beam  formed  in  SABF  is  compared  with  a  threshold 
value  provided  as  input  to  this  program  to  determine  when  the  threshold  is  first 
exceeded.  This  information  is  required  by  other  programs. 

E.6.2.3  Output 

There  are  two  forms  of  output  generated  by  SABF:  magnetic  tape  and 
printer. 

E. 6. 2. 3.1  Magnetic  Tape 

A  tape  consisting  of  pertinent  execution  parameters  and  all  beam  and 
filtered  beam  values  generated  is  produced  for  use  by  other  programs.  This 
tape  also  includes  the  value  of  the  center  seismometer  in  each  subarray  for 
each  sample  period  considered  in  the  processing.  The  logical  format  of  the 
tape  is  presented  on  a  sample  time  basis.  All  beams,  filtered  beams  and 
seismometer  values  formed  or  selected  per  sample  period  are  recorded  in  one 
logical  record.  The  total  number  of  logical  records  is  equal  to  the  rate  at 
which  the  beams  are  formed. 

E.6.2.3. 2  Printer 

Hard-copy  results  of  this  program  include  such  pertinent  run  parameters 
as  event  identification,  seismometers  considered,  the  number  of  samples  pro¬ 
cessed,  the  number  of  beams  formed,  the  time  delays  used  in  beamforming, 
filter  parameters  used,  and  the  points  at  which  the  beam  and  filtered  beam 
values  exceeded  the  threshold.  An  option  is  also  available  to  print  the  seis¬ 
mometer  and  beam  values  obtained  throughout  the  process. 

E.6.3  Program  Interaction 

This  program  requires  a  LASA  edited  tape  generated  by  the  LASA  Tape 
Edit  program.  The  filter  coefficients  used  are  obtained  from  the  results  of  the 
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Filter  programs.  Finally,  the  time  delays  used  are  provided  by  the  Threshold 
Time  Delay  program.  The  magnetic  tape  output  of  this  program  serves  as  input 
to  the  LASA  Beamformer  (LBF)  program  and  also  to  the  SABF  Power  Analysis 
program.  The  points  at  which  the  thresholds  are  exceeded  are  also  used  by  both 
of  these  programs. 


E.6.4  Program  Restrictions 

Data  restrictions  imposed  on  this  program  are  as  follows: 

Maximum  number  of  seismometers  per  subarray  -  25 

Maximum  number  of  beams  formed  per  subarray  -  6 

Highest  filter  degree -  6 

Maximum  time  history  of  seismometers  (sample  periods)  —  21 


The  above  data  restrictions  can  be  circumvented  for  individual  computer  runs 
if  trade -offs  with  other  data  restrictions  prove  effective. 

The  program  requires  that  the  number  of  seismometers  present  on  the 
input  tape  be  in  the  same  sequence  as  they  appear  on  the  raw  data  tapes. 

E.6.5  Comments 

This  program  is  written  in  FORTRAN  IV  except  for  a  subroutine  which 
reads  the  LASA  Edit  tape.  This  read  routine  is  written  in  MAP.  The  program 
currently  runs  on  the  IBM  7090  under  control  of  the  IBSYS  monitor.  Possible 
changes  to  this  program  include  the  following: 

a.  Provide  capability  to  vary  the  sampling  period 

b.  Provide  a  filter  on  the  seismometer  data  prior  to  their  use  in  the 
beamformer. 

c.  Provide  ability  to  process  subsets  of  the  seismometers  present  in 
each  subarray 

d.  Provide  ability  to  eliminate  any  individual  seismometer  in  the 
processing. 
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E.7  LASA  BEAMFORMER  (LBF) 


E.7.1  Purpose 

This  program  forms  a  central  LASA  beam  and  its  associated  neighboring 
beams.  The  beams  formed  by  this  program  will  then  be  available  for  further 
processing. 

E.7. 2  Description 

E.7 .2.1  Input 

The  output  tape  generated  by  the  Subarray  Beamformer  program  (SBF)  con¬ 
taining  both  unfiltered  and  filtered  beam  values  is  the  main  input  to  this  program. 
In  addition  to  this  beam  tape,  the  time  delays  required  to  form  the  central  beam 
and  the  neighboring  beams  for  both  unfiltered  and  filtered  beams  are  supplied 
by  card  input. 

E.7. 2.2  Computation 

The  beamforming  process  performed  by  this  program  is  the  same  as  the 
SABF  beamforming  process  except  that  instead  of  the  time  delayed  sum  of 
seismometer  data,  the  LASA  beams  are  formed  by  a  time  delayed  sum  of  sub¬ 
array  beams.  This  process  is  performed  first  for  the  central  beam  and  its 
neighboring  beams  using  unfiltered  subarray  beams  as  the  inputs,  and  then  for 
the  central  and  neighboring  beams  using  filtered  subarray  beams  as  the  inputs. 

Additional  processing  provides  the  sample  period  of  time  in  which  each  of  the 
LASA  beams  formed  exceeds  a  prescribed  threshold. 

E.7 .2. 3  Output 

This  program  produces  a  tape  of  LASA  beams  formed  over  a  selected  period 
of  time.  The  logical  format  of  the  tape  is  presented  on  a  sample  time  basis. 

All  beam  and  filtered  beams  formed  per  sample  period  are  recorded  in  one  logical 
record. 
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E.7.3  Program  Interaction 

The  subarray  beam  input  to  this  program  will  be  supplied  by  the  output  tape 
generated  by  the  SABF  program.  The  time  delays  for  the  neighboring  beams 
will  be  supplied  by  the  Neighboring  Beams  Time  Delay  Calculations  program. 

The  output  of  this  program,  the  LASA  beam  tape,  will  serve  as  input  to  the 
LBF  Power  Analysis  program  and  will  also  be  used  in  the  LASA  Beam  Pattern 
Display  process. 

E.7.4  Program  Restrictions 

At  the  present  time,  the  total  number  of  beams  this  program  can  form  is 
restricted  to  160.  The  maximum  time  history  for  each  beam  is  325  sample 
periods. 

E.7.5  Comments 

The  program  has  the  option  to  form  either  unfiltered  or  filtered  LASA  beams 
or  both.  The  type  of  LASA  beam,  filtered  or  unfiltered,  is  determined  by  the 
type  of  subarray  beam  used  in  the  beamforming  process  (filtered  or  unfiltered). 
No  filtering  is  done  in  the  LBF  program. 

The  program  is  written  in  FORTRAN  IV  for  the  IBM  7090. 

E.8  NEIGHBORING  BEAM  TIME  DELAY  CALCULATIONS 

E.8.1  Purpose 

This  program  generates  a  set  of  time  delays  to  form  a  close-packed  system 
of  beams  in  relation  to  a  particular  central  beam  being  investigated. 

E.8. 2  Description 

E.8 .2.1  Input 

The  main  input  to  this  program  is  the  location  of  the  source  of  the  event 
undergoing  analysis  which  is  approximated  by  the  location  in  inverse  velocity 
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space  of  the  centerline  of  the  best  detection  beam,  and  the  arrival  time  of  the 
wavefront  of  this  event  at  a  set  of  known  locations  (e.g.,  the  center  seismometer 
in  each  subarray  of  LASA). 


E.8.2.2  Computation 

A  set  of  neighboring  beams  related  to  the  centerline  of  the  detection  beam 
is  defined  by  causing  (Ux,  U  )  to  increment  so  as  to  form  neighboring  beams, 
and  limiting  the  neighborhood  by  prescribing  that 

(U  -  U  )2  +  (U  -  U  )2  <  R2 

I  '  x  xc7  '  y  yc7  I  — 

where: 


U  )  =  the  centerline  location  of  the  detection  beam  in 
yc7 

inverse-velocity  space 

R  =  the  radius  of  the  circle  of  the  selected  area  of 
detection  beam  uncertainly. 


A  closed-packed  system  may  be  defined  by  incrementing  (U  , 
the  rules: 


Uy)  according  to 


Ux  =uxc  +  ri3(m+0.5n) 
U  =  U  +  1.5  nr 

y  yc 


where  r  is  the  3.0  db  radius  of  an  event  beam  in  (Ux,  U^)  space,  and  m  and  n 
are  integers. 


til 

The  delay  at  the  k  subarray  whose  center  location  is  taken  as  (X^,  Y^,)  and 
for  a  beam  aimed  at  the  location  (Ux>  U  )  can  be  expressed  as  follows: 

t.  =  t.  -  (U  -  U  )  X,  +  (U  -  U  )  Y 
k  kc  x  xc7  k  y  yc  k 


where: 


*kc  =  -<*k  Uxc  +  Yk  V 

For  each  neighboring  beam  within  the  3.0  db  radius  of  event  circle,  (Ux, 
the  corresponding  time  delays  are  recorded. 


U  )  and 

y 
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E.8.3  Program  Interaction 

The  location  in  inverse  velocity  space  of  the  centerline  of  the  best  detection 

beam,  (U  ,  U  ),  is  obtained  from  the  output  of  the  Least  Squares  Wavefront 
xc  y  c 

program.  The  output  of  this  program,  the  sets  of  time  delays  and  centerline 
location  in  inverse -velocity  space  for  neighboring  beams  within  the  3.0  db  radius 
of  event  circle,  are  used  by  the  LASA  Beamformer  (LBF)  program. 

E.8.4  Program  Restrictions 

The  program  has  been  arbitrarily  restricted  to  21  subarrays. 

E.8.5  Comments 

The  program  is  written  in  FORTRAN  IV  for  the  IBM  System/360. 

E  .9  POWER  ANA  LYSIS  PROGRAMS 

E.9.1  Seismometer  Power  Analysis  Program 

E.9.1.1  Purpose 

This  program  calculates  both  noise  and  signal  power  at  the  seismometer 
level  for  any  particular  event  investigated. 

E.9.1. 2  Description 

E.9.1. 2.1  Input 

The  main  input  to  the  program  is  a  LASA  edited  tape  which  includes  all  the 
seismometer  channels  which  are  to  be  analyzed.  Additional  inputs  include  param¬ 
eters  specifying  the  starting  point  and  duration  for  both  the  noise  and  signal  power 
analysis,  and  the  relative  event  arrival  times  at  each  of  the  seismometers. 
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E. 9. 1.2. 2  Computation 


The  following  equation  represents  the  computation  to  determine  the  power 
in  the  seismometer  data  channels: 


where: 

N  =  number  of  samples  to  be  included  in  the  power  sum 
a  =  scale  factor 

ii  iL 

f.(n)  =  the  value  of  the  i  seismometer  at  the  n  sample  period 
P.  =  the  power  of  the  i1  seismometer. 

In  addition  to  the  power  calculations  on  an  individual  seismometer  basis, 
the  average  power  on  both  a  subarray  and  LASA  level  is  also  computed.  The 
variance  and  standard  deviation  of  the  seismometer  data  channels  are  also 
computed. 

E. 9. 1.2. 3  Outputs 

The  output  of  this  program  is  a  list  of  the  power  sums  (P.)  and  10  log  (P.) 
for  both  noise  and  signal.  The  average  power  values  described  previously  are 
also  recorded  along  with  the  pertinent  variance  and  standard  deviation  informa¬ 
tion. 


E.9.1.3  Program  Interaction 

Input  to  this  program  must  be  a  tape  recorded  in  the  LASA  Tape  Edit  format. 
System  gains  can  be  determined  by  comparing  the  results  of  this  program  with  the 
results  of  the  Subarray  and  LASA  beam  power  analysis  programs. 

E.9.1.4  Program  Restrictions 

The  noise  and  signal  intervals  must  be  separated  by  one  or  more  sample 
periods. 
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E.9.1.5  Comments 


The  program  is  written  in  FORTRAN  IV  for  the  IBM  7090. 

E.9.2  SABF  Power  Program 
E.9.2.1  Purpose 

This  program  calculates  noise  and  signal  power  present  in  the  beams  formed 
in  the  Subarray  Beamformer  program  (SABF).  This  data  may  then  be  compared 
with  a  similar  power  analysis  at  a  seismometer  level  and  LASA  beam  level  to 
indicate  overall  system  gain. 

E.9.2. 2  Description 

E.9.2. 2.1  Input 

The  main  input  to  this  program  is  the  tape  generated  by  the  SABF  program. 
This  tape  includes  unfiltered  and  filtered  beam  values  formed  for  some  specified 
period  of  time.  Other  input  values  include  parameters  indicating  points  at  which 
noise  and  signal  analysis  is  to  begin,  duration  of  both  noise  and  signal  power 
analysis,  and  the  points  at  which  each  beam  present  on  the  tape  exceeded  some 
prescribed  threshold  in  the  SABF  program. 

E.9.2. 2.2  Computation 

The  main  computations  in  this  program  are  the  power  calculations  and 
standard  deviation  performed  on  both  noise  and  signal  portions  of  the  filtered 
and  unfiltered  beams  present  on  the  input  tape. 

The  following  equation  represents  the  computation  to  determine  the  power 
of  a  beam: 
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where: 


N  =  number  of  samples  to  be  included  in  the  power  sum 
a  =  scale  factor 

th  th 

B.(n)  =  the  value  of  the  j  subarray  beam  at  the  n  sample 

J  th 

Pj  =  power  of  the  j  subarray  beam. 


In  addition  to  the  power  calculations  on  an  individual  beam  basis,  the  aver¬ 
age  beam  power  for  all  subarrays  present,  P,  is  also  calculated  and  the  associ¬ 
ated  standard  deviations.  This  process  is  done  first  for  the  unfiltered  and  then 
the  filtered  beams. 


E.9.2.2.3  Output 

The  output  of  this  program  is  a  list  of  the  individual  unfiltered  beam  power 
sums  (Pj)  and  10  log  (P^)  for  both  noise  and  signal.  In  addition,  the  average 
noise  and  signal  power  for  the  unfiltered  beams  with  the  associated  standard 
deviation  is  also  presented. 

The  same  output  is  recorded  for  the  filtered  beams. 

An  optional  output  is  provided  which  presents  the  intermediate  noise 
power  every  20  seconds  and  the  intermediate  signal  power  every  0.2 
seconds. 


E.9.2.3  Program  Interaction 

This  program  requires  an  SABF  output  tape  containing  subarray  beam 
values  over  some  specified  period  of  time.  In  addition,  the  output  of  the  SABF 
also  provides  data  indicating  the  start  of  the  event  for  each  beam  recorded  on 
the  SABF  output  tape. 

The  resultant  powers  calculated  in  this  program  can  be  compared  to  a 
similar  power  analysis  on  both  a  seismometer  level  and  a  IASA  beam  level  to 
indicate  overall  system  gain. 
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E.9.2.4  Program  Restrictions 

At  the  present  time,  the  program  is  designed  to  process  only  one  unfiltered 
and  one  filtered  beam  per  subarray.  The  program  also  expects  the  beam  input 
values  to  be  contained  on  a  single  tape.  The  selection  of  the  starting  points  for 
both  the  noise  and  signal  analysis,  as  well  as  the  duration  of  both,  must  be  done 
with  some  caution  to  ensure  that  the  required  amount  of  data  is  present  on  the 
beam  input  tape. 

E.9.2.5  Comments 

This  program  is  written  in  FORTRAN  IV  for  the  IBM  7090. 

E.10  THRESHOLD  TIME  DELAY 
E.10.1  Purpose 

This  program  provides  the  approximate  time  delays  required  to  form  the 
beams  in  the  Subarray  Beamformer  program.  No  attempt  is  made  to  seek 
maximum  gain  or  to  minimize  losses. 

E.10.2  Description 

E.  10. 2.1  Input 

The  main  input  to  this  program  is  a  selected  LASA  edited  tape  containing 
an  event  of  sufficient  magnitude  to  allow  for  the  successful  execution  of  the 
threshold  logic.  Other  pertinent  run  parameters,  such  as  the  threshold  value, 
are  provided  by  input  data  cards. 

E.10.2. 2.  Computation 

The  initial  computation  in  the  program  is  the  identification  of  the  sample 
periods  in  which  each  seismometer  present  on  the  input  tape  first  exceeds  a 
prescribed  threshold.  Once  each  seismometer  has  exceeded  the  threshold, 
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time  delays  (on  a  subarray  basis)  are  obtained  by  computing  the  difference  in 
sample  periods  between  the  time  each  seismometer  within  a  subarray  exceeded 
the  threshold  and  the  time  the  center  seismometer  in  the  subarray  exceeded  the 
threshold.  The  time  delay  for  the  i^  seismometer  in  the  subarray  is: 

TD(i,j)  =  THC(j)  -  TH(i,j)  +  T(j) 
where: 

THC(j) 

TH(i,  j) 

T  (j) 

E.10.2.3  Output 

The  sample  period  times  at  which  each  seismometer  exceeds  the  threshold 
as  well  as  the  time  delays  for  each  seismometer  relative  to  its  associated  sub¬ 
array  are  printed.  Additional  run  parameters  are  also  recorded. 

E.10.3  Program  Interaction 

The  time  delays  generated  by  this  program  are  used  in  the  Subarray  Beam- 
former  program  (SABF)  to  form  the  subarray  beams. 

E.10.4  Program  Restrictions 
None. 

E.10.5  Comments 

The  program  is  written  in  FORTRAN  IV  for  the  IBM  7090. 

Because  of  the  occurrence  of  corrupt  data,  the  thresholds  may  trigger 
thereby  generating  incorrect  time  delays. 

This  has  proved  useful  in  the  localization  of  the  "glitches"  and  the  time 
delays  are  easily  adjusted  by  investigation  of  the  data. 


=  the  sample  period  time  in  which  the  center  seismometer  in  the 
th 

j  subarray  exceeded  the  threshold 

th 

=  the  sample  period  time  in  which  the  1  seismometer  in  the 
j  n  subarray  exceeded  the  threshold 
=  sample  period  bias  for  the  j  subarray. 
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E.ll  SUBARRAY  PLOT 


E.ll.l  Purpose 

This  program  plots  seismometer  data  from  a  selected  subarray  for  a  vari¬ 
able  period  of  time.  The  resultant  plot  provides  a  visual  representation  of  the 
seismometer  activity  during  the  selected  period  of  time. 

E.11.2  Description 

The  Subarray  Plot  program  produces  a  magnetic  tape  for  a  CALCOMP  digi¬ 
tal  incremental  plotter.  The  plotter  produces  a  graphic  representation  of  the 
values  of  seismometers  within  a  specific  subarray  over  a  selected  period  of 
time. 

E. 11.2.1  Input 

Seismometer  data  is  obtained  from  a  LASA  edited  tape  while  control  param¬ 
eters  and  other  pertinent  data  such  as  the  selected  subarray  identification  is 
read  from  cards. 

E.ll. 2. 2  Computation 

The  program  reads  the  LASA  edited  tape  for  a  specified  period  of  time  and 
extracts  after  each  read  the  seismometer  values  of  the  selected  subarray. 

The  array  of  seismometer  data  obtained  over  the  prescribed  period  of  time  is 
then  placed  on  tape  by  the  CALCOMP  subroutines  for  subsequent  plotting. 

E.  11.2.3  Output 

Data  derived  from  the  input  tape  header  and  input  cards  are  printed  along 
with  execution  control  parameters  in  a  descriptive  format  to  provide  pertinent 
run  documentation.  The  magnetic  tapes  produced  by  the  program  are  in  the 
prescribed  format  for  off-line  plotting  on  the  CALCOMP  plotter. 


E.11.3  Program  Interaction 

The  program  requires  as  input  a  tape  generated  by  the  LASA  Tape  Edit 
program. 

E.11.4  Program  Restrictions 

Seismometer  values  are  read  and  plotted  in  22-second  blocks  with  a  limit 
of  6  blocks  (132  seconds)  per  plot.  The  capability  to  continue  the  plot  is  pro¬ 
vided  by  additional  passes  through  the  program.  The  maximum  number  of 
channels  (seismometers)  per  subarray  to  be  plotted  is  25. 

E.11.5  Comments 

The  program  is  written  in  FORTRAN  IV  for  the  IBM  7090.  Plotting  is  done 
on  a  CALCOMP  Model  563  plotter. 

E .  12  DISPLAY  PROGRAMS 

E.12.1  Display  Test  Tape  Program 

E.  12. 1.1  Purpose 

This  program  generates  four  tapes  to  check  the  Experimental  Display. 

The  display  is  described  in  Appendix  C. 

E.12.1. 2  Description 

E.12.1. 2.1  Input 

Input  to  the  program  consists  only  of  an  indication  of  which  tape  format  is 
to  be  produced.  A  fifth  possible  entry  requests  the  program  to  list  a  generated 
tape  in  a  format  which  approximates  that  of  the  display  face. 
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E.  12. 1.2. 2  Computation 

Based  on  the  input,  the  program  internally  sets  up  and  writes  out  the  desired 
tape  or  listing.  The  tapes  available  display  a  complete  face  of  a  single  intensity 
or  a  bar  pattern  consisting  of  all  16  possible  intensities.  The  single  intensity 
is  incremented  through  eight  possible  values  and  the  bar  pattern  is  rotated 
vertically  to  show  all  16  intensities  on  each  line  of  the  face.  Both  display  pat¬ 
terns  are  available  in  both  a  60-kc  and  a  30-kc  display  rate. 

E.  12. 1.2.3  Output 

The  program  outputs  are  tapes  and  listings  as  requested  by  the  operator. 

E.  12. 1.3  Program  Interaction 

This  program  stands  alone  and  is  independent  of  other  LASA  programs. 

E. 12. 1.4  Program  Restrictions 

Only  two  formats  can  be  produced.  Only  two  display  rates  (60  kc  and 
30  kc)  are  available. 

E. 12. 1.5  Comments 

This  program  was  implemented  in  the  Assembler  Language  of  the  Basic 
Programming  Support  for  the  IBM  System/360.  It  requires  a  machine  with  one 
tape  drive,  a  card  reader,  a  printer  and  a  console  typewriter.  The  last  require¬ 
ment  could  be  removed  by  using  parameter  cards  instead  of  operator  entries. 

E.12.2  LASA  Beam  Pattern  Display  Program 

E.  12 .2.1  Purpose 

The  LASA  Beam  Pattern  Display  program  prepares  input  tapes  for  use  on 
the  Experimental  Display  (see  Appendix  C). 
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E.12.2.2  Description 


E. 12.2. 2.1  Input 

Input  to  this  program  is  a  magnetic  tape  consisting  of  LASA  beam  sums. 

The  format  of  the  tape  is  defined  for  ease  of  generation  by  FORTRAN  and  ease 
of  acceptance  by  this  program.  The  tape  may  be  either  seven  or  nine  track. 
Other  program  inputs  are  parameter  cards  to  specify  the  action  desired  and 
address  cards  which  place  the  beam  sums  on  the  display  face.  This  technique 
allows  input  data  to  come  in  any  well  defined  order  as  long  as  the  address  cards 
are  constructed  in  the  same  order. 

E.12.2.2. 2  Computation 

The  beam  sums  are  read  and  converted  to  display  intensities  through  the 
use  of  a  table  look-up  procedure.  This  allows  a  look  at  any  range  of  the  beam 
signals  in  detail  while  suppressing  all  other  portions. 

E.12.2.2 .3  Output 

Two  outputs  are  generated  by  the  program.  First  is  a  display  tape  which 
contains  from  1  to  99  copies  of  a  display  sequence.  The  number  of  copies  is  an 
input  parameter.  The  second  output  is  a  listing  of  selected  records  of  the  tape 
in  a  format  similar  to  that  displayed  on  the  display  device. 

E.12.2.3  Program  Interaction 

The  magnetic  tape  input  comes  from  the  LASA  Beamformer  program  after 
a  reformatting  process.  The  address  cards  are  produced  by  the  Display  Address 
Card  program. 

E.12.2.4  Program  Restrictions 

Only  two  display  duration  times  (50  -  100  milliseconds)  are  available.  Only 
one  output  tape  can  be  produced  although  the  input  may  extend  over  nine  reels  of 
tape,  if  necessary. 
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E.12.2.5  Comments 


This  program  was  written  in  Assembler  Language  for  the  IBM  System /360. 

It  is  designed  to  operate  using  the  Basic  Programming  Support  System.  Input 
requests  may  be  stacked  and  tape  or  listing  requests  may  be  intermixed  in  any 
order. 

E.12.3  LASA  Display  Address  Card  Program 
E. 12. 3.1  Purpose 

This  program  generates  address  cards  for  use  as  input  to  the  LASA  Beam 
Pattern  Display  program. 

E. 12. 3.2  Description 

E.12.3. 2.1  Input 

The  input  to  this  program  consists  of  the  deck  of  cards  generated  by  the 
neighboring  beam  program.  This  deck  contains  the  coordinates  of  each  beam 
in  velocity  space. 

E.  12.3.2.2  Computation 

Each  beam  is  placed  in  an  internal  display  matrix  using  the  inverse  velocity 
space  coordinates.  The  first  beam  in  sequence  is  placed  in  the  center  and  all 
other  beams  are  placed  relative  to  it.  The  method  employed  for  placing  the  beams 
creates  holes  in  the  pattern  such  that  beams  are  not  displayed  in  contiguous 
display  positions.  A  compaction  process  is  then  used  with  the  location  of  the 
center  beam  remaining  constant. 

E.12.3 .2. 3  Output 

A  punched  card  deck  of  addresses  and  a  listing  of  velocity  coordinates  is 
produced.  The  punched  cards  are  suitable  for  input  to  the  LASA  Beam  Pattern 
Display  program. 
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E. 12.3.4  Program  Interaction 

The  input  of  this  program  is  the  output  of  the  neighboring  beams  program. 
The  output  of  this  program  serves  as  part  of  the  input  to  the  LASA  Beam  Pattern 
Display  program. 

E. 12. 3.5  Program  Restrictions 

Any  velocity  coordinates  which  initially  place  the  beam  off  the  display  face 
are  ignored  and  treated  as  zero. 

E.  12.3.6  Comments 

The  program  was  written  in  FORTRAN  IV  for  use  on  an  IBM  System/360. 
The  program  requires  33,544  bytes  of  storage. 

E.13  CONVERTED  STUDY  PROGRAMS 
E.13.1  Seismic  Ray-Tracing  Program 

E.13. 1.1  Purpose 

This  program  calculates  the  angular  range,  travel  time,  and  inverse  phase 
speed  for  both  the  direct  and  surface  reflected  rays,  and  repeats  this  calculation 
for  arbitrary  equal  increments  of  Snell's  constant  over  any  desired  range  of 
values. 

E.13. 1.2  Description 
E.13. 1.2.1  Input 

This  program  requires  a  number  of  layers,  N,  and  for  each  layer  the  radii 
and  velocities  at  upper  and  lower  boundaries.  The  source  position  is  defined  to 
be  at  the  upper  boundary  of  any  layer  with  the  receiver  at  the  surface. 
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E.  13. 1.2. 2  Computation 

This  program  was  written  under  the  LASA  Signal  Processing  Study*  and  has 
been  converted  to  System/360  FORTRAN  IV. 

E.  13. 1.2. 3  Output 

For  each  layer,  the  velocity  gradient  and  zero-radius  velocity  are  printed 
along  with  travel  time.  Also,  angular  distance  of  both  direct  and  reflected  ray 
inverse  phase  speed,  and  radius  of  deepest  penetration  for  each  value  of  Snell’s 
constant  is  recorded. 

E.  13.1.3  Program  Interaction 

Relation  of  inverse  phase  speed  to  range  is  used  in  the  Inverse  Velocity 
Space  Mapping  program. 

E. 13. 1.4  Program  Restrictions 

In  its  current  form,  the  program  is  restricted  to  95  layers.  The  source 
position  must  be  at  the  upper  boundary  of  one  of  the  95  layers. 

E. 13. 1.5  Comments 

The  program  is  written  in  FORTRAN  IV  for  the  IBM  System  /360. 

E.13.2  Least  Squares  Wavefront  with  Range  Correction  Program 
E.  13. 2.1  Purpose 

This  program  calculates  the  direction  and  speed,  and  thereby  the  approxi¬ 
mate  azimuth  and  range,  of  the  wavefront  corresponding  to  a  given  set  of  arrival 
times  at  a  set  of  known  locations  by  calculating  the  best-fitting  plane  wavefront 
corresponding  to  these  arrivals.  An  attempt  is  made  to  correct  the  arrival 
times  for  curvature  as  a  function  of  range. 

*IBM  Final  Report,  "LASA  Signal  Processing  Study,"  Contract  No.  SD-296, 

July  15,  1965. 
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E.  13.2 .2  Description 


E.  13 .2. 2.1  Input 

A  set  of  arrival  times  at  a  set  of  known  locations  is  required  as  input. 

t 

E.  13 .2 .2 .2  Computation 

This  program  was  written  under  the  LASA  Signal  Processing  Study*  and  has 
|  been  converted  to  System/360  FORTRAN  IV. 

E. 13.2. 2. 3  Output 

The  program  provides  azimuth  and  range  to  event  and  the  deviation  in  seconds 
of  the  actual  arrivals  from  the  best-fitting  plane  wavefront. 

E.13.2.3  Program  Interactions 

The  punched  output  from  this  program  (deviations  of  the  actual  arrivals  from 
the  best -fitting  plane  wavefront)  is  used  as  input  to  the  Seismic  Steering  Delay 
1  Anomalies  program. 

E.13.2.4  Program  Restrictions 

,  To  determine  a  plane  wavefront,  the  arrival  times  at  three  different  locations 

must  be  available  for  the  event.  The  program  is  arbitrarily  restricted  to  a 
maximum  of  21  subarrays. 

E.13.2.5  Comments 

This  program  was  written  in  FORTRAN  IV  for  the  IBM  System/360. 


♦IBM  Final  Report,  "LASA  Signal  Processing  Study,"  Contract  No.  SD-296, 
July  15,  1965. 
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E.13.3  Seismic  Steering  Delay  Anomalies  Program 
E.  13. 3.1  Purpose 

This  program  calculates  the  average  and  standard  deviation  of  the  deviation 
from  the  best-fitting  plane  wavefront  at  each  seismometer*  from  an  average  for 
that  seismometer,  and  calculates  the  loss  in  db  suffered  by  each  wavefront 
because  of  mis -steering. 

E.13.3. 2  Description 

E.13.3. 2.1  Input 

The  program  accepts  deviations  in  seconds  of  the  actual  arrivals  from  the 
best -fitting  plane  wavefront  for  each  wavefront  and  each  seismometer. 

E.13.3. 2. 2  Computation 

The  average  and  standard  deviation  of  wave  arrival  deviations  from  plane 
are  calculated  for  each  seismometer,  along  with  the  loss  in  db  due  to  mis- 
steering  for  each  wavefront.  The  loss  is  found  from 

Loss  =  (10  log10  e)  (27rf  a)2 

where: 

a  =  the  standard  deviation 
f  =  the  frequency. 

E.13.3. 2. 3  Output 

The  average  arrival  time  deviations,  differences  of  deviations  from  average, 
and  db  loss  are  recorded. 


*  Seismometer,  as  used  here,  may  also  be  interpreted  as  the  center  of  a 
subarray. 
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E.13.3.3  Program  Interaction 


This  program  accepts  as  card  input  the  deviations  of  the  actual  arrivals 


from  the  best-fitting  plane  wavefront  from  the  Least  Squares  Wavefront  program. 


E.13.3.4  Program  Restrictions 

This  program  arbitrarily  limits  the  number  of  subarrays  to  21  and  limits 
the  number  of  wavefronts  to  200.  There  must  be  at  least  three  time  deviations 
for  each  wavefront. 

E.13.3.5  Comments 

This  program  was  written  in  FORTRAN  IV  for  the  IBM  System /360. 

E.14  CROSS  COVARIANCE  OF  SEISMIC  DATA  CHANNELS 
E.14.1  Purpose 

The  program  computes  and  plots  the  cross -covariance  of  two  seismic  data 
input  channels  (seismometers  or  beams). 

E.14. 2  Description 

E.  14.2.1  Input 

The  program  accepts  digitized  data  in  several  formats.  Control  cards 
specify  the  seismometers  and  segments  of  information  to  be  processed. 

E.14. 2. 2  Computation 

The  cross-covariance  formulae,  where  A.,  i  =  1,...,  n,  and  B^,  i  =  1,...,  m, 
are  two  sets  of  data,  and  n  <  m  are: 


T 


A.  •  B.+j  _1  j  /T  -(Average  A  •  Average  B);  j  =  1,...,  m 


i=l 
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where: 


T  =  minimum  (n,  m  -  j  +  1), 


Average  A 


Ai  /n; 


Average  B 


and 


S. 

J 


A.  j  •  B.  J  /T  -(Average  A  -  Average  B);  j  =  1 


n 


where: 


T  =  n  -j  +  1. 


(The  above  is  a  modification  of  the  method  of  J.  R.  Healy  and  B.  P.  Bogert  as 
found  in  the  "System/360  Scientific  Subroutine  Package  Programmer's  Manual," 
IBM  Corporation,  Form  H20-0205. ) 

E.14.2.3  Output 

The  output  consists  of  a  CALCOMP  plotter  tape  which  results  in  an  annotated 
graph  of  each  seismometer  and  the  cross -covariance  curve. 

E.14.3  Program  Interaction 
None. 

E.14.4  Program  Restrictions 

The  selected  segment  of  any  seismometer  cannot  be  greater  than  2400  observa¬ 
tions.  (This  would  be  two  minutes  of  LASA  data  if  every  sample  period  is  considered.) 

E.14.5  Comments 

The  program  is  written  in  IBM  7090  FORTRAN  IV  but  could  be  converted 
for  use  on  another  machine  if  CALCOMP  plotter  subroutines  were  available  for 
that  computer. 
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E.15  POWER  SPECTRA 


E.15.1  Purpose 

The  Power  Spectra  program  calculates  the  Fourier  transform  of  seismom¬ 
eter  traces  each  consisting  of  N  real  data  samples  obtained  from  a  LASA  edit 
tape  or  tapes.  The  number  of  samples  is  taken  to  be  a  power  of  2.  Fourier 
transforms  are  obtained  for  two  seismometer  traces  at  a  time.  The  program 
will  also  operate  on  a  single  seismometer  trace  but  with  some  loss  in  efficiency. 

E.15. 2  Description 

E.  15.2. 1  Input 

Input  for  the  program  consists  of:  (a)  an  800  BPI  FORTRAN  compatible 
LASA  edit  tape  containing  the  seismometer  traces,  and  (b)  at  least  one  control 
card  specifying  two  seismometer  channel  numbers,  m  for  the  sample  size 
N  =  2m  to  be  taken  from  the  trace  of  each  seismometer  specified,  and  the  number 
of  time  samples  to  be  skipped  before  processing.  The  second  of  the  two 
seismometer  channel  numbers  is  optional. 

E.15. 2. 2  Computation 

Given  a  set  of  N  =  2m  real  data  samples  X(j),  j  =  •  *,  N  - 1,  the  program 

computes  the  set  of  N  complex  coefficients  A(k),  k  =  0,  1,  •  •  •  *,  N  -  1  of  the  skewed 
complex  Fourier  series : 


N-l 


k=0 


where: 


W 


i 
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The  Aik)  coefficients  are  found  by  the  inverse  formula 


N-l 

A(k)=^-  ^  X(j)  W_jk  k  =  0,  1  N-l 

j=0 

Tie  Cooley-Tukey  algorithm11  is  employed  in  the  calculations 

t 

E.  15.2.3  Output 

Output  of  the  program  consists  of  an  off-line  printout  of:  (a)  the  time  and 
channel  data  appearing  in  the  header  record  of  the  LASA  edit  tape,  (b)  tabulation 
of  the  data  samples  taken  from  each  seismometer  trace,  and  (c)  a  tabulation  of 
the  complex  Fourier  series  coefficients  A(k),  k  =  0,  1,  •  •  •,  N/2.  Because  the 
data  samples  are  real,  the  remaining  A(k)  coefficients  are  not  printed  for 
A(N-k)  =  A(k)  k=  1,  2,  •  •  •,  -g-  -1,  where  A(k)  means  the  conjugate  of  A(k). 

E .  15 .3  Program  Interaction 
Note. 

E.15.4  Program  Restrictions 

The  sample  size  taken  from  a  given  seismometer  trace  must  be  a  power  of 

13 

2  not  exceeding  8192  =2  or  less  than  2. 

E .  15 .5  Comments 

The  program  is  written  in  FORTRAN  IV  and  MAP  for  the  IBM  7090. 


*  J.  W.  Cooley  and  J.  W.  Tukey,  "An  Algorithm  for  the  Machine  Calculation  of 
Complex  Fourier  Series,"  Mathematics  of  Computation,  Vol.  19,  April  1965,  p  297. 
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Appendix  F 


COMMENTS  AND  RECOMMENDATIONS  ON  THE  IBM  LASA  STUDY 

F.l  INTRODUCTION 

'jTiis  appendix*  presents  comments  on  the  following  topics  related  to  the 
IBM  LASA  study  and  system  efforts: 

a  .  Developing  the  Steering  Delays  (Section  F.2) 
t .  Consideration  of  Post-P  Arrivals  (Section  F  .3) 

c.  Processing  and  Processor  Configuration  (Section  F.4) 

d.  Array  Synthesis  (Section  F. 5) 

e.  Source  Location  and  Use  of  Several  LASAs  (Section  F.6) 

f .  Problems  for  Further  Research  (Section  F.7) 

g.  LASA  Seismic  Handbook  (Section  F.8). 

F.2  DEVELOPING  THE  STEERING  DELAYS 

Specification  of  a  first  set  of  beam  steering  delays  to  be  programmed  into  the 
detection  processor  as  soon  as  it  becomes  available  is  perhaps  the  most  urgent 
task  to  be  undertaken.  It  is  clear  that  such  delays  must  at  first  be  based  essen¬ 
tially  on  experimental  data  of  seismic  wave  fronts  actually  received  at  LASA. 

The  current  IBM  effort  to  relate  the  delays  to  those  corresponding  to  idealized 
plane  wave  fronts  appears  to  be  a  reasonable  approach  to  this  problem.  This 
approach  has  a  possible  advantage  over  the  standard  method  using  Jeffreys-Bullen 
travel  times  with  local  station  corrections  in  that  it  separates  the  problem  of 
obtaining  the  steering  delays  from  the  problem  of  source  location.  For  each 
wave  iront,  the  steering  delays  are  determined  in  such  a  way  as  to  minimize  the 
residuals  from  an  idealized  front  rather  than  to  conform  to  a  separately  computed 
source  location  which  in  any  case  will  often  become  available  only  some  time 
after  the  event. 

*This  ippendix  is  based  on  notes  by  A.  F.  Gangi,  Associate  Professor,  Department 
of  Geology  and  Geophysics,  MIT,  and  Consultant  to  IBM  LASA  Project. 
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The  present  plane  wave  front  approach  uses  a  least  squares  plane  wave  front 

to  classify  the  steering  delays  according  to  inverse  velocity  components  (W  ,  W  ) 

x  y 

and  lists  the  deviations  from  that  plane.  There  might  be  an  advantage  to  using 
a  best  fitting  quadratic  wave  front  instead,  so  as  to  permit  study  of  the  two  wave 
curvatures.  If  these  curvatures,  in  addition  to  the  velocities,  should  be  found  to 
behave  in  a  reasonably  predictable  manner,  then  the  quadratic  reference  wave 
fronts  would  of  course  permit  smaller  residuals,  and  hence,  improved  steering 
than  the  plane  wave  method. 

Initially,  all  steering  delays  must  be  based  on  an  entirely  empirical  formula¬ 
tion,  using  experimental  delays.  Based  on  plane  (or  quadratic)  wave  fronts,  the 
delays  must  be  adjusted  for  some  average  value  of  the  residual  at  each  subarray. 
However,  it  is  very  important  for  several  reasons  to  also  begin  to  base  this 
information  on  a  more  exact  theoretical  foundation  in  terms  of  the  earth's  struc¬ 
ture  of  propagation  velocities  which  in  turn  defines  travel  times  and  phase  veloc¬ 
ities.  This  is  what  in  effect  has  already  been  done  when  the  Jeffreys-Bullen  Tables 
were  constructed.  However,  with  the  advent  of  large  array  technology,  it  should 
become  possible  to  reach  into  the  velocity  profile  with  greater  detail  and  precision, 
based  on  wave  velocity  and  curvature  measurements  as  well  as  travel  times . 

The  most  immediate  reason  why  this  theoretical  basis  must  be  established 
is  to  provide  a  good  method  for  steering  to  those  parts  of  the  earth  from  which 
seismic  events  of  adequate  strength  are  relatively  scarce.  Another  reason  is  that 
an  accurate  knowledge  of  the  earth's  velocity  structure  will  permit  more  accurate 
utilization  of  later  arrivals  to  identify  events  and  pinpoint  their  locations  and  times 
of  origm.  Finally,  a  better  knowledge  of  the  earth's  structure  is  of  itself  of  great 
scientific  interest  and  can  in  the  longer  run  lead  to  a  more  systematic  understanding 
of  geophysical  questions . 

F.3  CONSIDERATION  OF  POST-P  ARRIVALS 

Se:.smic  arrivals  following  the  initial  P  wave  will  be  received,  and  must, 
therefore,  be  identified  even  if  their  further  analysis  will  not  always  be  required. 

In  addiion,  these  later  arrivals,  when  grouped  and  identified,  will  aid  in  the 
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discrimination  problem  and  may  improve  event  location  accuracy.  Arrivals  of 
importance  will  be  pP,  PcP,  S,  and  others.  pP  and  to  some  extent  PcP  should 
be  particularly  useful  in  source  depth  estimation.  Amplitude  of  the  various  phases 
should  be  compared  to  begin  to  accumulate  information  on  seismic  source  radia¬ 
tion  patterns  in  P  and  S  modes.  Such  patterns  may  become  useful  discriminants. 

F.4  PROCESSING  AND  PROCESSOR  CONFIGURATION 

The  necessity  for  a  "post-detection  or  pre-event"  processing  function  is 
disci.ssed,  although  this  function  might  be  included  in  the  event  processor.  The 
problem  is  to  weed  out  in  a  logical  manner  the  many  arrivals  the  detection 
processor  will  detect,  including  post-P  and  side-lobe  detections.  The  following 
funct.ons  in  particular  should  be  considered. 

a.  Detect  and  identify  pP  to  obtain  an  estimate  of  source  depth,  if  any. 

b.  Detect  and  identify  PcP  to  obtain  further  information  on  source 
depth  and  distance  to  the  source. 

c .  Detect  and  identify  other  arrivals  with  relatively  large  amplitudes . 

d.  Decide  which  events  and  which  arrivals  should  be  subjected  to 
further  processing  by  the  event  processor . 

Cnee  some  of  the  arrivals  have  been  identified,  beams  should  be  formed 
at  correct  times  and  velocities  to  definitely  identify  as  many  of  the  original 
detections  as  practicable  and,  where  necessary,  to  produce  clean  wave  forms 
and  olher  required  event  processor  outputs. 

If  it  turns  out  that  event  beamforming  and  event  beam  selections  place  an 
excessive  load  on  beamforming  capability,  it  may  be  desirable  to  form  successively 
closei  together  event  beams,  always  dividing  up  the  remaining  area  into  a  few 
beams  and  selecting  the  best  one,  rather  than  immediately  covering  the  whole 
area  with  many  closely  spread  beams . 

A  suggestion  is  made  to  investigate  an  adaptive  technique  using  correlation 
to  sharpen  an  event  beam  by  adaptively  adjusting  its  steering  delays. 

Tie  Seismic  Bulletin  to  be  published  should  include  not  only  the  classical 
seism.c  station  outputs,  but  also  the  noise  statistics,  the  signal  phase  speeds, 
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and  perhaps  the  wave  front  curvatures  as  measured  by  the  array. 

The  possibility  of  improving  signal-to-noise  ratio  by  vectorially  combining 
the  three  outputs  of  the  three-component  seismometer  should  also  be  investigated. 

F.5  ARRAY  SYNTHESIS 

It  may  be  worthwhile  to  make  an  array  geometry  synthesis  program  based 
on  the  superposition  principle.  The  array  patterns  of  a  triangle,  square,  pentagon, 
and  hexagon  should  be  made.  Then,  beam  patterns  of  arrays  composed  of  combina¬ 
tions  of  these  simple  geometries  can  be  readily  calculated  by  taking  into  account 
the  scale  change  on  a  pattern  in  wave  number  space  of  a  change  in  the  array  size 
in  physical  space.  Also,  a  rotation  in  physical  space  corresponds  to  a  rotation 
in  wave  number  space.  This  should  be  a  great  help  in  the  synthesis  of  patterns. 

It  is  not  necessary  to  limit  oneself  to  the  regular  polygons  in  performing  this 
study,  and  it  also  is  not  necessary  to  synthesize  array  patterns  using  just  one 
type  of  polygon  for  the  whole  array. 

F.6  SOURCE  LOCATION  AND  USE  OF  SEVERAL  LASAS 

The  use  of  more  than  one  LASA  can  be  expected  to  somewhat  improve  loca¬ 
tion  accuracy  over  a  single  LASA  because : 

a.  One  of  several  LASAs  might  be  in  a  more  advantageous  location 
relative  to  the  event. 

b.  An  improvement  of  the  order  of  ^  N  when  using  N  LASAs  should 
follow  statistically. 

c .  If  four  or  more  LASAs  are  used,  then  world-wide  base  lines  for 
calculating  source  location  become  available,  greatly  improving 
the  location  accuracy.  However,  the  same  effect  can  be  achieved 
using  single  seismometers  in  place  of  LASAs  except  that  the  signal- 
to-noise  ratio  would  be  lower,  and  hence,  the  detection  identifica¬ 
tion  and  localization  threshold  would  necessarily  be  at  a  higher 
signal  level. 
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F.7  PROBLEMS  FOR  FURTHER  RESEARCH 


So  me  problems  which  should  be  investigated  by  using  the  LASA  are : 

a.  The  velocity  distribution  in  the  earth 

b.  The  attenuation  of  waves  in  the  earth 

c.  Earthquake  statistics 

d.  Variation  of  earthquake  spectra  with  magnitude 

e .  Radiation  patterns  of  earthquakes 

f.  Microseismic  noise  statistics 

g.  Weather  or  storm  location 

h.  Local  geological  structure 

i.  Earthquake  prediction. 

A  better  knowledge  of  the  velocity  distribution  in  the  earth  can  be  obtained 
from  the  measured  arrival  times,  horizontal  velocities,  and  curvatures  of  the 
waves  incident  on  the  array.  This  problem  is  closely  related  to  that  of  deter¬ 
mining  the  proper  steering  delays  to  form  beams.  In  addition,  the  amplitudes 
of  the  wives  can  also  provide  information  regarding  the  density  distribution  in 
the  earth.  The  variation  of  the  velocity  and  density  distributions  with  azimuth 
would  delineate  the  lateral  inhomogeneities  in  the  properties  of  the  earth's 
interior 

From  the  difference  in  the  spectra  of  the  direct  P  arrival  and  later  arrivals, 
some  measure  of  the  attenuation  in  the  earth  can  be  obtained.  Before  these 
measurements  can  be  made  with  accuracy,  it  will  be  necessary  to  determine  the 
effects  of  the  local  geology  and  the  source  radiation  patterns  on  the  spectra  of  the 
received  waves . 

With  the  automatic  detection,  localization,  depth  determination,  and  (to  some 
degree)  magnitude  determination  ability  of  LASA,  along  with  its  sensitivity  to 
small  magnitude  earthquakes,  it  should  be  very  easy  to  obtain  very  good  earth¬ 
quake  statistics  in  a  short  time  and  without  the  normal  drudgery  involved  in 
obtaining  these  statistics.  These  data  are  also  of  value  to  the  operational  system 
to  determine  some  of  the  operating  parameters. 

Because  the  signals  received  at  LASA  will  be  digitized  and  available  on  mag¬ 
netic  tape,  it  will  be  reasonably  simple  to  obtain  the  spectra  of  these  signals. 
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The  correlation  of  the  characteristics  of  these  spectra  with  other  earthquake 
parameters  would  lead  to  a  better  understanding  of  the  earthquake  source  mech¬ 
anism. 

With  a  world-wide  distribution  of  LASAs  or  even  with  the  detection  of  later 
arrivals  from  a  single  LASA,  it  should  be  possible  to  obtain  information  on  earth¬ 
quake  radiation  patterns  on  a  routine  basis .  This  too  would  lead  to  a  better 
understanding  of  the  source  mechanism  of  the  underlying  causes  of  tectonic 
activity  in  the  earth. 

The  study  of  microseismic  noise  statistics  would  lead  to  a  better  knowledge 
of  the  local  geology  and  its  effects  on  the  propagation  of  microseismic  noise.  In 
addition,  because  microseismic  noise  is  believed  to  be  caused  by  weather  (storms) 
at  sea  or  near  the  coasts,  this  study  of  the  microseisms  could  be  used  to  locate 
and  track  such  storms.  The  information  obtained  here  would  also  be  of  value  to 
the  operational  system. 

The  local  geological  structure  can  be  studied  by  a  knowledge  of  the  station 
corrections  at  each  site  and  their  variations  with  distance  and  azimuth.  Further 
information  on  the  local  or  region  structure  can  be  obtained  by  looking  at  small, 
local  events  such  as  local  quakes  or  quarry  blasts.  In  this  sense,  the  array  can 
be  used  as  a  tool  to  perform  refraction  profiles . 

Finally,  it  is  possible  that  LASA  can  be  used  as  a  tool  for  earthquake  predic¬ 
tion.  This  study  would  be  performed  in  conjunction  with  the  earthquake  statistics 
studies  because  it  has  been  indicated  that  earthquake  activity  (small  events) 
increase  prior  to  larger  events.  This  theory  is  not  very  well  substantiated  and 
LASA  would  be  a  convenient  tool  to  verify  or  refute  the  hypothesis . 

F.8  LASA  SEISMIC  HANDBOOK 

Aaother  suggestion  is  to  begin  work  on  a  LASA  Seismic  Handbook  in  which 
characteristics  of  the  array  are  kept.  This  would  include  such  things  as  beam 
pattens,  processing  gain,  signal  loss  due  to  timing  errors,  noise  statistics  and 
their  temporal  (seasonal)  variation.  Earthquake  statistics,  maps  of  seismic 
regions  as  determined  by  LASA  and  many  other  kinds  of  data.  It  would  not  be 
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necessary  to  explain  the  given  results  or  to  explain  how  they  were  obtained.  It 
would  suff  ice  to  give  references  to  the  original  papers  and  reports  from  which 
the  data  were  taken.  In  this  way,  all  the  pertinent  data  of  value  in  the  operation 
of  LASA  would  be  at  one's  fingertips  in  condensed  form. 
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